This code generates the figures from the outputs of the model.
binder, where is the .yml file for the conda virtual environment.
data, where are the outputs of the model NEMO-LIM3, the outputs of the model MITgcm, and the observations from the Qikiqtarjuaq ice camps 2015 and 2016.
figures_progress, where are the figures.
This notebook was run in a conda virtual environment. The following commands will automatically create the environment and launch the Jupyter notebook:
cd <my_path>/wintertime/
conda env create -f binder/environment.yml
conda activate wintertime
cd utils/python/MITgcmutils
python setup.py bdist_egg
cd build
pip install MITgcmutils
cd ../../../..
jupyter notebook
Once the environment wintertime is created, the Jupyter notebook can be launched simply with:
cd <my_path>/wintertime/
conda activate wintertime
jupyter notebook
import cartopy
import cartopy.crs as ccrs
import collections
import datetime
import gsw
import gudinfo
import IPython.display
import netcdf_tools
import math
import matplotlib as mpl
import matplotlib.colors
import matplotlib.image
import matplotlib.patches
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.transforms
import numpy as np
import os
import os.path
import pandas as pd
import PIL
import PIL.Image
import read_mitgcm
import scipy
import scipy.interpolate
from scipy.io import loadmat
import vstats
plt.close("all")
Before running the script, download the file containing the IBCAO bathymetry in the current directory: https://bit.ly/3vqJAHt (134 MB).
'''latlon2xy convenience function to convert latlon to chosen projection,
and figure object'''
def latlon2xy(lon,lat,proj):
xyz = proj.transform_points(ccrs.PlateCarree(),lon,lat).squeeze()
x,y,_ = np.split(xyz,3,axis=len(xyz.shape)-1)
x[np.isinf(x)]=np.nan
y[np.isinf(y)]=np.nan
return x.squeeze(),y.squeeze()
# STATIONS
Station=collections.namedtuple(
'Station','lon lat marker markeredgecolor markerfacecolor markersize')
station1=Station(-60.3919999999999,68.3096,'o','none','red',5)
station2=Station(-61.463,68.0358,'o','none','red',5)
station3=Station(-62.355,67.8668,'o','none','red',5)
station4=Station(-63.5639999999999,67.54,'o','none','red',5)
station5=Station(-63.783,67.482,'*','black','yellow',10) # Qikiqtarjuaq SIC
station6=Station(-64.638,67.239,'o','none','red',5)
stations=[station1,station2,station3,station4,station5,station6]
plotbathymetry=True
with plt.style.context('mplstyles/maps.mplstyle'):
fig=plt.figure(figsize=(10, 16))
# --- CCGS AMUNDSEN 2018, LEG 2B
proj_am=ccrs.LambertConformal(central_longitude=-63,
central_latitude=67.5,
standard_parallels=(67,68))
ax_am=fig.add_subplot(223,projection = proj_am)
ax_am.coastlines(resolution='10m')
ax_am.add_feature(cartopy.feature.LAND)
if not plotbathymetry:
ax_am.add_feature(cartopy.feature.OCEAN)
MapExtent = [-66,-60,66,69]
ax_am.set_extent(MapExtent, ccrs.PlateCarree())
gl = ax_am.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
linewidth=1, color='gray', alpha=0.5,
linestyle='--')
gl.xlocator = mticker.FixedLocator(np.arange(-68,-58))
gl.ylocator = mticker.FixedLocator(np.arange(65,70))
# stations
for station in stations:
xstation,ystation = latlon2xy(np.array([station.lon]),
np.array([station.lat]),
proj_am)
marker=station.marker
markerfacecolor=station.markerfacecolor
markersize=station.markersize
markeredgecolor=station.markeredgecolor
ax_am.plot(xstation,ystation,marker,markeredgecolor=markeredgecolor,
markerfacecolor=markerfacecolor,markersize=markersize)
# longitude labels
y0,_ = ax_am.get_ylim()
for lon in [-66,-64,-62,-60]:
x0,_ = proj_am.transform_point(lon,MapExtent[2]-1,
src_crs = ccrs.PlateCarree())
plt.text(x0,y0-2e4,'{:2d}$^\circ$W'.format(-lon),
horizontalalignment = 'center',verticalalignment='top')
# latitude labels
x0,_ = ax_am.get_xlim()
someLons = np.arange(-80,-50)
for lat in np.arange(66,69+1):
# interpolate latitude circle to map boundary
xyzArray = proj_am.transform_points(ccrs.PlateCarree(),
someLons,lat*np.ones_like(someLons))
x = xyzArray[:,0]
y = xyzArray[:,1]
y0 = np.interp(x0,x,y)
plt.text(x0-2e4,y0,'{:2d}$^\circ$N'.format(lat),
horizontalalignment = 'right',verticalalignment='center')
if plotbathymetry:
# color bar for bathymetry
cmap = plt.cm.get_cmap('Greys')
cmap = cmap.from_list('Custom cmap',
[cmap(i) for i in range(150,50,-1)], 100)
cmap.set_over('white')
bounds = [0,500,1000,2000]
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
# IBCAO bathymetry
if plotbathymetry:
top = loadmat('data/DataS1_observations_IBCAO_1min_bathy.mat')
lon_ibcao,lat_ibcao=np.meshgrid(top['x'],top['y'])
z_ibcao=-top['IBCAO_1min']
hb_am=ax_am.contourf(*latlon2xy(lon_ibcao,lat_ibcao,proj_am),
z_ibcao,levels=[0,500,1000,2000],
cmap=cmap,norm=norm,extend='max')
axcb_am= fig.colorbar(hb_am,ax=ax_am,orientation='horizontal')
axcb_am.ax.set_title('Water depth (m)')
axcb_am.ax.set_position([0.70,-0.10, 0.25,0.20])
# labels for stations
for i,station in enumerate(stations):
xstation,ystation = latlon2xy(np.array([station.lon]),
np.array([station.lat]),
proj_am)
ax_am.text(xstation+3e3,ystation+1e4,i+1,
horizontalalignment = 'center',
verticalalignment='center')
# text
b_lon,b_lat = -66.5, 69.1,
ax_am.text(b_lon,b_lat,'(b)',
horizontalalignment = 'left',fontsize=30,
transform=ccrs.PlateCarree())
# --- BAFFIN BAY
proj_bb=ccrs.LambertConformal(central_longitude=-60,
central_latitude=70,
standard_parallels=(65,75))
ax_bb=fig.add_subplot(221,projection=proj_bb)
ax_bb.coastlines(resolution='50m')
ax_bb.add_feature(cartopy.feature.LAND)
if not plotbathymetry:
ax_bb.add_feature(cartopy.feature.OCEAN)
MapExtent = [-76,-48,62,79]
ax_bb.set_extent(MapExtent, ccrs.PlateCarree())
gl = ax_bb.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
linewidth=1, color='gray', alpha=0.5,
linestyle='--')
gl.xlocator = mticker.FixedLocator([-80,-75,-70,-65,-60,
-55,-50,-45,-40,-35])
gl.ylocator = mticker.FixedLocator(np.arange(58,84,3))
# stations
for station in stations:
xstation,ystation = latlon2xy(np.array([station.lon]),
np.array([station.lat]),
proj_bb)
marker=station.marker
markerfacecolor=station.markerfacecolor
markersize=station.markersize
markeredgecolor=station.markeredgecolor
ax_bb.plot(xstation,ystation,marker,markeredgecolor=markeredgecolor,
markerfacecolor=markerfacecolor,markersize=markersize)
# longitude labels
y0,_ = ax_bb.get_ylim()
for lon in [-70,-60,-50]:
x0,_ = proj_bb.transform_point(lon,MapExtent[2]-1,
src_crs = ccrs.PlateCarree())
plt.text(x0,y0-2e4,'{:2d}$^\circ$W'.format(-lon),
horizontalalignment = 'center',verticalalignment='top',
fontsize=20)
# latitude labels
x0,_ = ax_bb.get_xlim()
someLons = np.arange(-120,-40,1)
for lat in np.arange(64,78,3):
# interpolate latitude circle to map boundary
xyzArray = proj_bb.transform_points(ccrs.PlateCarree(),
someLons,lat*np.ones_like(someLons))
x = xyzArray[:,0]
y = xyzArray[:,1]
y0 = np.interp(x0,x,y)
plt.text(x0-2e4,y0,'{:2d}$^\circ$N'.format(lat),
horizontalalignment = 'right',verticalalignment='center')
# IBCAO bathymetry
if plotbathymetry:
hb_bb=ax_bb.contourf(*latlon2xy(lon_ibcao,lat_ibcao,proj_bb),
z_ibcao,levels=[0,500,1000,2000],
cmap=cmap,norm=norm,extend='max')
# text
a_lon,a_lat = -98, 77.5
ax_bb.text(a_lon,a_lat,'(a)',
horizontalalignment = 'left',fontsize=30,
transform=ccrs.PlateCarree())
greenland_lon,greenland_lat = -38, 77
ax_bb.text(greenland_lon,greenland_lat,'Greenland',
horizontalalignment = 'right',
transform=ccrs.PlateCarree())
baffinbay_lon,baffinbay_lat = -60, 72
ax_bb.text(baffinbay_lon,baffinbay_lat,'Baffin\nBay',
horizontalalignment = 'right',
transform=ccrs.PlateCarree())
# rectangle
extentStudyarea = ax_am.get_extent()
vertices = np.array([
[extentStudyarea[0], extentStudyarea[2]], # left, bottom
[extentStudyarea[0], extentStudyarea[3]], # left, top
[extentStudyarea[1], extentStudyarea[3]], # right, top
[extentStudyarea[1], extentStudyarea[2]], # right, bottom
[extentStudyarea[0], extentStudyarea[2]]
])
xVert,yVert = vertices.transpose()
xBox = np.interp(np.arange(0,4,.02),np.arange(0,4.1,1),xVert)
yBox = np.interp(np.arange(0,4,.02),np.arange(0,4.1,1),yVert)
xyzBoxLatlon = ccrs.PlateCarree().transform_points(ax_am.projection,
xBox,yBox)
xBoxLatlon = xyzBoxLatlon[:,0]
yBoxLatlon = xyzBoxLatlon[:,1]
xyzBoxTransformed = proj_bb.transform_points(ccrs.PlateCarree(),
xBoxLatlon,yBoxLatlon)
xBoxTransformed = xyzBoxTransformed[:,0]
yBoxTransformed = xyzBoxTransformed[:,1]
ax_bb.plot(xBoxTransformed,yBoxTransformed,'r')
ax_bb.relim()
# --- WORLD MAP
proj_wm = ccrs.NearsidePerspective(central_longitude=-60,
central_latitude=60)
ax_wm = fig.add_subplot(222,projection = proj_wm)
ax_wm.set_global()
ax_wm.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')
ax_wm.add_feature(cartopy.feature.OCEAN)
gl = ax_wm.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
linewidth=1, color='gray', alpha=0.5,
linestyle='--')
extentStudyarea = ax_bb.get_extent()
vertices = np.array([
[extentStudyarea[0], extentStudyarea[2]], # left, bottom
[extentStudyarea[0], extentStudyarea[3]], # left, top
[extentStudyarea[1], extentStudyarea[3]], # right, top
[extentStudyarea[1], extentStudyarea[2]], # right, bottom
[extentStudyarea[0], extentStudyarea[2]]
])
xVert,yVert = vertices.transpose()
xBox = np.interp(np.arange(0,4,.02),np.arange(0,4.1,1),xVert)
yBox = np.interp(np.arange(0,4,.02),np.arange(0,4.1,1),yVert)
xyzBoxLatlon = ccrs.PlateCarree().transform_points(ax_bb.projection,
xBox,yBox)
xBoxLatlon = xyzBoxLatlon[:,0]
yBoxLatlon = xyzBoxLatlon[:,1]
xyzBoxTransformed = proj_wm.transform_points(ccrs.PlateCarree(),
xBoxLatlon,yBoxLatlon)
xBoxTransformed = xyzBoxTransformed[:,0]
yBoxTransformed = xyzBoxTransformed[:,1]
ax_wm.plot(xBoxTransformed,yBoxTransformed,'r')
ax_wm.set_position( [0.60,0.60 ,0.38,0.38])
ax_bb.set_position( [0.10,0.53 ,0.60,0.42])
ax_am.set_position( [0.10,0.05 ,0.60,0.42])
# --- SAVE
if plotbathymetry:
plt.savefig('figures_progress/maps.png')
else:
plt.savefig('figures_progress/maps_nobathy.png')
years=range(1,11)
nbyears=len(years)
iyeartempo=nbyears
first_year=np.arange(0,365)
last_year=np.arange(365*(nbyears-1)+0,365*(nbyears-1)+365)
array1d_iT1y_iT=last_year
first_year366=np.arange(0,366) # first year for heatmaps
drF is the r cell face separation, meaning the thickness of each depth layer (in m)
It corresponds to delR on https://mitgcm.readthedocs.io/en/latest/getting_started/getting_started.html#grid
gridfile='data/DataS8_output_mitgcm/exp0/grid.t001.nc';
drF=netcdf_tools.read_netcdf(gridfile,'drF')
RC is the r coordiante of cell center (in m). We switch sign to have positive depths.
RC=-netcdf_tools.read_netcdf(gridfile,'RC')
RF is the r coordinate of cell interface (in m). We switch sign to have positive depths.
RF=-netcdf_tools.read_netcdf(gridfile,'RF')
RF_above151=np.array(RF[RF<151])
Break-up was on 18 July 2016 (Oziel et al., 2019 in Elem. Sci. Anth.).
obsice=np.empty(365)
obsice[:]=np.NaN
obsice[0:199]=1
obsice[199]=0
Day of in situ sea ice break-up:
obs_breakup=next(i for i,v in enumerate(obsice) if v < 0.50)
The file Ice thickness from https://www.seanoe.org/data/00487/59892/ contains $in\ situ$ snow and ice thickness at the Qikiqtarjuaq sea ice camps 2015 and 2016 ($67.48^\circ N$, $-63.79^\circ E$). Variable 'sample_thickness_cm_average' is the snow or ice thickness (cm).
def load_thickness():
fname='data/DataS2_observations_Qikiqtarjuaq/66407.csv'
df = (
pd.read_csv(
fname,
dtype={
'mission': 'string',
'date': str,
'date_time': str,
'latitude': np.float32,
'longitude': np.float32,
'sample_type': 'category',
'sample_thickness_cm_average': np.float32,
'sample_thickness_cm_sd': np.float32,
'pi': 'category'
},
parse_dates=['date','date_time']
)
)
df['doy']=df['date'].apply(lambda x:x.timetuple().tm_yday)
df['sample_thickness_m_average']=df['sample_thickness_cm_average'] \
.apply(lambda x:x/100)
return df
thickness_df=load_thickness()
We select the observations of snow thickness at the ice camp 2016.
snow_thickness_df=thickness_df[(thickness_df.mission=='ice_camp_2016') \
& (thickness_df.sample_type=='snow')]
We select the observations of ice thickness at the ice camp 2016.
ice_thickness_df=thickness_df[(thickness_df.mission=='ice_camp_2016') \
& (thickness_df.sample_type=='ice')]
The file C-OPS from https://www.seanoe.org/data/00487/59892/ contains $in\ situ$ underwater PAR at the Qikiqtarjuaq sea ice camps 2015 and 2016 ($67.48^\circ N$, $-63.79^\circ E$). Variable 'averaged_par_d_fit_daily_ein_m_2_day_1' is the downwelling PAR (mol photons $\mathrm{ m^{-2} }$ $\mathrm{ d^{-1} }$).
def load_par():
fname='data/DataS2_observations_Qikiqtarjuaq/66399.csv'
df = (
pd.read_csv(
fname,
dtype={
'mission': 'string',
'station': 'category',
'date': str,
'snow_thickness': 'category',
'depth_m': np.float32,
'pi': 'category',
'averaged_par_d_fit_muein_m_2_s_1': np.float32,
'averaged_par_d_fit_percent_percent': np.float32,
'averaged_par_d_fit_daily_ein_m_2_day_1': np.float32,
'averaged_par_d_noon1hloc_ein_m_2_1h_1': np.float32,
'averaged_par_d_p1h_ein_m_2_1h_1': np.float32,
'averaged_par_d_p3h_ein_m_2_3h_1': np.float32,
'averaged_par_d_p24h_ein_m_2_day_1': np.float32,
'averaged_par_d_p48h_ein_m_2_48h_1': np.float32
},
parse_dates=['date']
)
.rename(
columns=dict(
depth_m='depth',
averaged_par_d_fit_daily_ein_m_2_day_1='par',
)
)
)
df['doy']=df['date'].apply(lambda x:x.timetuple().tm_yday)
return df
par_df=load_par()
par_df=par_df[par_df.mission=='ice_camp_2016']
find negative PAR
par_df['par'][par_df['par']<0]
Series([], Name: par, dtype: float32)
find duplicates
g=par_df.groupby('doy')['depth'].value_counts()
g.where(g>1).dropna()
Series([], Name: depth, dtype: float64)
The first value of simulated PAR had a depth=0 m even when there was sea ice. The first value of observed PAR had a depth equals to the thickness of sea ice in the water. For example, if there was 1 m of sea ice in the water, the first value of observed PAR would be at depth=1 m. We changed that for a depth=0 m to be coherent with the simulated values. Indeed, we wanted to compare observed and simulated PAR at the same distance from the bottom of sea ice and not at the same distance from the water surface.
parwater_df=(par_df.dropna(subset=['par'])).copy()
for doy in pd.unique(parwater_df['doy']):
ice_thickness=parwater_df.loc[(parwater_df.doy==doy),('depth')].min() # m
parwater_df.loc[parwater_df.doy==doy,'depth']\
=parwater_df.loc[parwater_df.doy==doy,'depth']-ice_thickness
array1d_iT_obsPAR3m=np.empty(365)
array1d_iT_obsPAR3m[:]=np.NaN
for doy in pd.unique(parwater_df['doy']):
y_interp=scipy.interpolate.interp1d(
parwater_df[parwater_df.doy==doy]['depth'],
parwater_df[parwater_df.doy==doy]['par'])
array1d_iT_obsPAR3m[doy-1]=y_interp(3)
array1d_iT_obsisolume=np.empty(365)
array1d_iT_obsisolume[:]=np.NaN
c_isolume=0.415 # mol photons m^-2 d^-1
for doy in pd.unique(parwater_df['doy']):
parwater_onedoy_df=parwater_df[parwater_df.doy==doy]
if parwater_onedoy_df.iloc[0,:]['par'] > c_isolume:
idx_isolume=vstats.find_idx_nearest \
(array=parwater_onedoy_df['par'],
value=c_isolume)
isolume=parwater_onedoy_df['depth'].iloc[idx_isolume]
array1d_iT_obsisolume[doy-1]=isolume
The file Nutrients from https://www.seanoe.org/data/00487/59892/ contains $in\ situ$ nutrients at the Qikiqtarjuaq sea ice camps 2015 and 2016 ($67.48^\circ N$, $-63.79^\circ E$). Variables 'no3_um_l', 'po4_um_l', and 'sioh4_um_l' are the nitrate, phosphate, and silicic acid concentrations, respectively ($\mathrm{ \mu }$M).
def load_nutrients():
fname='data/DataS2_observations_Qikiqtarjuaq/66412.csv'
df = (
pd.read_csv(
fname,
dtype={
'mission': 'string',
'date': str,
'date_time': str,
'latitude': np.float32,
'longitude': np.float32,
'sample_type': 'category',
'sample_source': 'category',
'snow_thickness': 'category',
'code_operation': 'category',
'station': 'category',
'station_type': 'category',
'cast': np.float32,
'bottle': np.float32,
'bottom_depth_m': np.float32,
'target_depth_m': 'string',
'depth_m': np.float32,
'filter_type': 'category',
'quality_flag': 'category',
'dilution_factor': np.float32,
'no3_um_l': np.float32,
'no2_um_l': np.float32,
'no2_and_no3': np.float32,
'po4_um_l': np.float32,
'sioh4_um_l': np.float32,
'n_ice_samples_ice_camp_2016': np.float32,
'method': 'category',
'pi': 'category'
},
parse_dates=['date','date_time']
)
.rename(
columns=dict(
depth_m='depth',
no3_um_l='no3',
po4_um_l='po4',
sioh4_um_l='sioh4'
)
)
)
df['doy']=df['date'].apply(lambda x:x.timetuple().tm_yday)
return df
nutrients_df=load_nutrients()
We select the observations of nitrate at the ice camp 2016 in the water column made with a GFF filter.
nutrients_df=nutrients_df[(nutrients_df.mission=='ice_camp_2016') \
& (nutrients_df.sample_type=='water') \
& (nutrients_df.filter_type=='gff')]
Find duplicates.
g=nutrients_df.groupby('doy')['depth'].value_counts()
g.where(g>1).dropna()
Series([], Name: depth, dtype: float64)
No duplicates.
array1d_iT_obsnitrate3m=np.empty(365)
array1d_iT_obsnitrate3m[:]=np.NaN
array1d_iT_obssilica3m=np.empty(365)
array1d_iT_obssilica3m[:]=np.NaN
for doy in pd.unique(nutrients_df['doy']):
y_interp=scipy.interpolate.interp1d(
nutrients_df[nutrients_df.doy==doy]['depth'],
nutrients_df[nutrients_df.doy==doy]['no3'])
array1d_iT_obsnitrate3m[doy-1]=y_interp(3)
for doy in pd.unique(nutrients_df['doy']):
y_interp=scipy.interpolate.interp1d(
nutrients_df[nutrients_df.doy==doy]['depth'],
nutrients_df[nutrients_df.doy==doy]['sioh4'])
array1d_iT_obssilica3m[doy-1]=y_interp(3)
in $\mathrm{ \mu M=mmol\ Fe\ m^{-3} }$
fname='data/DataS5_observations_davis_strait.csv'
fet_df=pd.read_csv(
fname,
dtype={
'doy':np.int64,
'depth':np.float32,
'Fe_D_CONC_BOTTLE_nmol_kg':np.float32
}
)
fet_df['depth']=-fet_df['depth']
# fet_df
# nmol kg^-1 ===> mmol m^-3
fet_df['Fe_D_CONC_BOTTLE_nmol_kg']=fet_df['Fe_D_CONC_BOTTLE_nmol_kg']/1000
fet_df.rename({'Fe_D_CONC_BOTTLE_nmol_kg':'Fe_D_CONC_BOTTLE_mmol_m3'},
axis='columns',
inplace=True)
# fet_df
The file Pigments, nutrients, Chlorophyll a and Phaeopigments (concentration) from https://www.seanoe.org/data/00487/59892/ contains $in\ situ$ chlorophyll a at the Qikiqtarjuaq sea ice camps 2015 and 2016 ($67.48^\circ N$, $-63.79^\circ E$). Variables 'conc_mg_m3' is the chlorophyll a concentration measured by HPLC (mg Chl $\mathrm{ m^{-3} }$).
def load_pigment():
fname='data/DataS2_observations_Qikiqtarjuaq/66417.csv'
df = (
pd.read_csv(
fname,
dtype={
'mission': 'string',
'date': str,
'date_time': str,
'latitude': np.float32,
'longitude': np.float32,
'sample_id': 'string',
'sample_code': 'string',
'sample_type': 'category',
'sample_source': 'category',
'snow_thickness': 'category',
'leg': 'category',
'code_operation': 'category',
'station': 'category',
'station_type': 'category',
'cast': np.float32,
'bottle': np.float32,
'target_depth_m': 'string',
'depth_m': np.float32,
'pigment': 'category',
'conc_mg_m2': np.float32,
'conc_ug_l': np.float32,
'conc_mg_m3': np.float32,
'filtered_vol_l': np.float32,
'vol_after_melt_l': np.float32,
'ice_added_vol_of_water': np.float32,
'ice_total_vol_water': np.float32,
'ice_calculated_vol_ice': np.float32,
'dilution_factor': np.float32,
'quality_control_peaks': 'category',
'nb_core': np.float32,
'nb_core_9cm_dia': np.float32,
'core_area_9cm_dia': np.float32,
'nb_core_14cm_dia': np.float32,
'core_area_14cm_dia': np.float32,
'area_m2': np.float32,
'description_of_sample': 'category',
'n': np.float32,
'observations': 'category',
'method': 'category',
'pi': 'category'
},
parse_dates=['date','date_time']
)
.rename(
columns=dict(
depth_m='depth',
conc_ug_l='chlfluo',
conc_mg_m3='chlHPLC'
)
)
)
df['doy']=df['date'].apply(lambda x:x.timetuple().tm_yday)
return df
chlHPLC_df=load_pigment()
We select the observations of chlorophyll a concentrations at the ice camp 2016 measured by HPLC.
chlHPLC_df=chlHPLC_df[(chlHPLC_df.mission=='ice_camp_2016') \
& (chlHPLC_df.sample_type=='water') \
& (chlHPLC_df.pigment=='Total Chlorophyll a') \
& (chlHPLC_df.method=='HPLC')]
remove negative chlorophyll a
chlHPLC_df['chlHPLC'][ chlHPLC_df['chlHPLC']<0 ]
Series([], Name: chlHPLC, dtype: float32)
find duplicates
g=chlHPLC_df.groupby('doy')['depth'].value_counts()
g.where(g>1).dropna()
doy depth
146 1.5 2.0
10.0 2.0
Name: depth, dtype: float64
inspect duplicates for doy=146
chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC','doy']][(chlHPLC_df.doy==146) | (chlHPLC_df.doy==148)]
chlHPLC_oneday_df
| depth | chlHPLC | doy | |
|---|---|---|---|
| 20357 | 0.5 | 0.1056 | 146 |
| 20358 | 1.5 | 0.0899 | 146 |
| 20359 | 5.0 | 0.0906 | 146 |
| 20360 | 10.0 | 0.0905 | 146 |
| 20361 | 20.0 | 0.1011 | 146 |
| 20362 | 40.0 | 0.0633 | 146 |
| 20363 | 0.5 | 0.0971 | 148 |
| 20364 | 1.5 | 0.1483 | 146 |
| 20365 | 5.0 | 0.1086 | 148 |
| 20366 | 10.0 | 0.0890 | 146 |
| 20367 | 20.0 | 0.0735 | 148 |
| 20368 | 40.0 | 0.0399 | 148 |
By human judgement, I replace the doy of rows 20364 and 20366 with 148.
chlHPLC_df.loc[[20364,20366],'doy']=148
chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC','doy']][(chlHPLC_df.doy==146) | (chlHPLC_df.doy==148)]
chlHPLC_oneday_df
| depth | chlHPLC | doy | |
|---|---|---|---|
| 20357 | 0.5 | 0.1056 | 146 |
| 20358 | 1.5 | 0.0899 | 146 |
| 20359 | 5.0 | 0.0906 | 146 |
| 20360 | 10.0 | 0.0905 | 146 |
| 20361 | 20.0 | 0.1011 | 146 |
| 20362 | 40.0 | 0.0633 | 146 |
| 20363 | 0.5 | 0.0971 | 148 |
| 20364 | 1.5 | 0.1483 | 148 |
| 20365 | 5.0 | 0.1086 | 148 |
| 20366 | 10.0 | 0.0890 | 148 |
| 20367 | 20.0 | 0.0735 | 148 |
| 20368 | 40.0 | 0.0399 | 148 |
Identify days of year for which the number of observations is less than or equal to 2. There are not measurements at enough depths to be able to calculate a meaningful vertical integration.
g=chlHPLC_df.groupby('doy')['doy'].value_counts()
g.where(g<=2).dropna()
doy doy 120 120 2.0 Name: doy, dtype: float64
chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC','doy']][(chlHPLC_df.doy==120)]
chlHPLC_oneday_df
| depth | chlHPLC | doy | |
|---|---|---|---|
| 20299 | 1.5 | 0.0355 | 120 |
| 20300 | 5.0 | 0.0253 | 120 |
By human judgement, I drop these observations.
chlHPLC_df=chlHPLC_df.drop(chlHPLC_df[(chlHPLC_df.doy==120)].index)
in $\mathrm{ mg\ Chl\ m^{-2} }$
#doy: 123
depths_case1 =np.array([ 1.5, 5, 10, 20, 40])
weights_case1=np.array([ 3.25,4.25, 7.5, 15, 10])
#doys: 125,127
depths_case2 =np.array([0.5, 1.5, 5, 20, 40])
weights_case2=np.array([1, 2.25,9.25, 17.5, 10])
#doys:130,132,134,137,139,144,146,148,151,153,155,158,160,162,165,167,169,172,176
depths_case3 =np.array([0.5,1.5, 5, 10, 20, 40])
weights_case3=np.array([1, 2.25, 4.25, 7.5, 15, 10])
#doy: 141
depths_case4 =np.array([0.5,1.5, 5, 10, 40])
weights_case4=np.array([1, 2.25, 4.25,17.5, 15])
#doy: 174
depths_case5 =np.array([0.5,1.5, 10, 20, 40])
weights_case5=np.array([1, 4.75, 9.25,15, 10])
#doys: 179,181,183,186,188,190
depths_case6 =np.array([0.5,1.5, 10, 20,30, 60])
weights_case6=np.array([1, 4.75, 9.25,10,15, 0])
#doys: 193,195,200
depths_case7 =np.array([ 1.5, 10, 20,30, 60, 75])
weights_case7=np.array([ 5.75, 9.25,10,15, 0, 0])
#doys: 197,204
depths_case8 =np.array([ 1.5, 10, 20,30, 45, 60])
weights_case8=np.array([ 5.75, 9.25,10,12.5, 2.5,0])
nT=365
array1d_iT_obsvint40mchlHPLC=np.full(nT,np.NaN)
for iT in np.arange(0,nT):
obsvintchlHPLC=np.NaN
chlHPLC=[]
doy=iT+1
chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC']][(chlHPLC_df.doy==doy)]
depths=chlHPLC_oneday_df['depth']
if set(depths_case1).issubset(set(depths)):
for depth in depths_case1:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case1)
elif set(depths_case2).issubset(set(depths)):
for depth in depths_case2:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case2)
elif set(depths_case3).issubset(set(depths)):
for depth in depths_case3:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case3)
elif set(depths_case4).issubset(set(depths)):
for depth in depths_case4:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case4)
elif set(depths_case5).issubset(set(depths)):
for depth in depths_case5:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case5)
elif set(depths_case6).issubset(set(depths)):
for depth in depths_case6:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case6)
elif set(depths_case7).issubset(set(depths)):
for depth in depths_case7:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case7)
elif set(depths_case8).issubset(set(depths)):
for depth in depths_case8:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case8)
array1d_iT_obsvint40mchlHPLC[iT]=obsvintchlHPLC
in $\mathrm{ mg\ Chla\ m^{-2} }$
#doy: 123
depths_case1 =np.array([ 1.5, 5, 10, 20, 40])
weights_case1=np.array([ 3.25,4.25, 7.5, 15, 70])
#doys: 125,127
depths_case2 =np.array([0.5, 1.5, 5, 20, 40])
weights_case2=np.array([1, 2.25,9.25, 17.5, 70])
#doys:130,132,134,137,139,144,146,148,151,153,155,158,160,162,165,167,169,172,176
depths_case3 =np.array([0.5,1.5, 5, 10, 20, 40])
weights_case3=np.array([1, 2.25, 4.25, 7.5, 15, 70])
#doy: 141
depths_case4 =np.array([0.5,1.5, 5, 10, 40])
weights_case4=np.array([1, 2.25, 4.25,17.5, 75])
#doy: 174
depths_case5 =np.array([0.5,1.5, 10, 20, 40])
weights_case5=np.array([1, 4.75, 9.25,15, 70])
#doys: 179,181,183,186,188,190
depths_case6 =np.array([0.5,1.5, 10, 20,30, 60])
weights_case6=np.array([1, 4.75, 9.25,10,20, 55])
#doys: 193,195,200
depths_case7 =np.array([ 1.5, 10, 20,30, 60, 75])
weights_case7=np.array([ 5.75, 9.25,10,20, 22.5,32.5])
#doys: 197,204
depths_case8 =np.array([ 1.5, 10, 20,30, 45,60])
weights_case8=np.array([ 5.75, 9.25,10,12.5,15,47.5])
nT=365
array1d_iT_obsvintchlHPLC=np.full(nT,np.NaN)
for iT in np.arange(0,nT):
obsvintchlHPLC=np.NaN
chlHPLC=[]
doy=iT+1
chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC']][(chlHPLC_df.doy==doy)]
depths=chlHPLC_oneday_df['depth']
if set(depths_case1).issubset(set(depths)):
for depth in depths_case1:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case1)
elif set(depths_case2).issubset(set(depths)):
for depth in depths_case2:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case2)
elif set(depths_case3).issubset(set(depths)):
for depth in depths_case3:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case3)
elif set(depths_case4).issubset(set(depths)):
for depth in depths_case4:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case4)
elif set(depths_case5).issubset(set(depths)):
for depth in depths_case5:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case5)
elif set(depths_case6).issubset(set(depths)):
for depth in depths_case6:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case6)
elif set(depths_case7).issubset(set(depths)):
for depth in depths_case7:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case7)
elif set(depths_case8).issubset(set(depths)):
for depth in depths_case8:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case8)
array1d_iT_obsvintchlHPLC[iT]=obsvintchlHPLC
obsphotofile='data/DataS7_literature.csv'
obsphoto_df=pd.read_csv(obsphotofile)
obsPCmax_df=obsphoto_df[obsphoto_df.parameter=='PCmax']
obsalphachl_df=obsphoto_df[obsphoto_df.parameter=='alphastar']
obsEk_df=obsphoto_df[obsphoto_df.parameter=='Ek']
The file siarea.nemo.2016.365.32bits.bin contains data from a simulation of NEMO3.6 LIM3.6 configured by Gaetan Olivier (UBO) for the Qikiqtarjuaq sea ice camp location ($67.48^\circ N$, $-63.79^\circ E$) in 2016. It contains sea ice concentration (unitless).
ice=np.fromfile('data/DataS4_output_nemo_lim3/siarea.nemo.2016.365.32bits.bin',
dtype='>f')
Day of simulated sea ice freeze-up:
sim_freezeup=next(i for i,v in reversed(list(enumerate(ice))) if v < 0.50)
The file Ice_d.nc is from a simulation of NEMO3.6 LIM3.6 configured by Gaetan Olivier (UBO) for the Qikiqtarjuaq sea ice camp location ($67.48^\circ N$, $-63.79^\circ E$) from 2013 to 2017. The variable snow_volume contains the snow volume (m). The variable ice_volume contains the ice volume (m). The variables are at indices (1,1) meaning centre of grid point was used.
Snow volume is the weighted average of snow thickness averaged by the sea ice concentration. For example, if there is 50% of sea ice concentration and 2 m of snow on the sea ice-covered part of the pixel, snow volume will be 1 m. The expression 'snow thickness' in the labels of the plots corresponds to the variable 'snow volume' in NEMO-LIM3 and the NetCDF file Ice_d.nc.
ice_file='data/DataS4_output_nemo_lim3/Ice_d.nc'
From 2013 to 2017.
array1d_iT_snow_volumefull=\
netcdf_tools.read_netcdf(ice_file,'snow volume')
array1d_iT_snow_volumefull=array1d_iT_snow_volumefull[:,1,1]
From 2016.
array1d_iT_snow_volume=array1d_iT_snow_volumefull[365*3:365*4]
Ice volume is the weighted average of ice thickness averaged by the sea ice concentration. For example, if there is 50% of sea ice concentration and 2 m of ice on the sea ice-covered part of the pixel, ice volume will be 1 m. The expression 'ice thickness' in the labels of the plots corresponds to the variable 'ice volume' in NEMO-LIM3 and the NetCDF file Ice_d.nc.
From 2013 to 2017.
array1d_iT_ice_volumefull=\
netcdf_tools.read_netcdf(ice_file,'ice volume')
array1d_iT_ice_volumefull=array1d_iT_ice_volumefull[:,1,1]
From 2016.
array1d_iT_ice_volume=array1d_iT_ice_volumefull[365*3:365*4]
indir='data/DataS8_output_mitgcm/exp0'
info=gudinfo.RunInfo(indir)
esd=info.f['esd'].to_numpy()
stod=86400 # seconds to day
pcmax=info.f['pcmax'].to_numpy() # s^-1
pcmax=pcmax*stod # d^-1
ksatno3=info.f['ksatno3'].to_numpy()
stod=86400 # seconds to day
gmax=info.f['grazemax'].to_numpy() # s^-1
gmax=gmax*stod # d^-1
The file car.0000000000.t001.nc was the output of the default (EXP-0) simulation.
The file car.0000000000.t001.nc (180 MB) was too large for GitHub. It is available on my Google Drive.
carfile='data/DataS8_output_mitgcm/exp0/car.0000000000.t001.nc'
($\mathrm{ mol\ photons\ m^{-2}\ d^{-1} }$)
stod=86400 # 3600*24 s per d
array2d_idepth_iT_PARfullumolEonm2s=netcdf_tools.read_netcdf(carfile, 'PAR').squeeze().transpose()
if np.all(array2d_idepth_iT_PARfullumolEonm2s==0):
#define GUD_PARUICE
array2d_idepth_iT_PAR_icefullumolEonm2s=netcdf_tools.read_netcdf(carfile, 'PAR_ice').squeeze().transpose()
array2d_idepth_iT_PAR_icefull=array2d_idepth_iT_PAR_icefullumolEonm2s*1e-6*stod
array2d_idepth_iT_PAR_owfullumolEonm2s=netcdf_tools.read_netcdf(carfile, 'PAR_ow').squeeze().transpose()
array2d_idepth_iT_PAR_owfull=array2d_idepth_iT_PAR_owfullumolEonm2s*1e-6*stod
(ndepths,ndays)=array2d_idepth_iT_PAR_icefull.shape
array2d_idepth_iT_icefull=np.tile(ice,(ndepths,ndays//365))
array2d_idepth_iT_meanPARfull=array2d_idepth_iT_PAR_owfull*(1-array2d_idepth_iT_icefull)\
+array2d_idepth_iT_PAR_icefull*array2d_idepth_iT_icefull
else:
#undef GUD_PARUICE
array2d_idepth_iT_meanPARfull=array2d_idepth_iT_PARfullumolEonm2s*1e-6*stod
array2d_idepth_iT_meanPAR=array2d_idepth_iT_meanPARfull[:,array1d_iT1y_iT]
(ndepths,ndays)=array2d_idepth_iT_meanPAR.shape
array1d_iT_modPAR3m=np.zeros(ndays)
x=3
xp=[RC[2],RC[3]]
assert np.all(np.diff(xp)>0), "xp is in decreasing order: %r %r" % (xp[0],xp[1])
for iT in first_year:
fp =[array2d_idepth_iT_meanPAR[3,iT],array2d_idepth_iT_meanPAR[2,iT]]
PAR3m=np.interp(x,xp,fp)
array1d_iT_modPAR3m[iT]=PAR3m
($\mathrm{ \mu mol\ photons\ m^{-2}\ s^{-1} }$)
array2d_idepth_iT_meanPARumolEonm2s=array2d_idepth_iT_meanPAR/1e-6/stod
c_isolume=0.415 # mol photons m^-2 d^-1
(ndepths,ndays)=array2d_idepth_iT_meanPAR.shape
array1d_iT_modisolume=np.empty(ndays)
array1d_iT_modisolume[:]=np.NaN
for iT in first_year:
array1d_idepth_meanPAR=array2d_idepth_iT_meanPAR[:,iT]
if array1d_idepth_meanPAR[0]>c_isolume:
idx_isolume=vstats.find_idx_nearest(array=array1d_idepth_meanPAR,
value=c_isolume)
isolume=RC[idx_isolume]
array1d_iT_modisolume[iT]=isolume
array2d_idepth_iT_modno3full=netcdf_tools.read_netcdf(carfile, 'TRAC04').squeeze().transpose()
array2d_idepth_iT_modno3=array2d_idepth_iT_modno3full[:,array1d_iT1y_iT]
array2d_idepth_iT_modsioh4full=netcdf_tools.read_netcdf(carfile, 'TRAC06').squeeze().transpose()
array2d_idepth_iT_modsioh4=array2d_idepth_iT_modsioh4full[:,array1d_iT1y_iT]
array2d_idepth_iT_modpo4full=netcdf_tools.read_netcdf(carfile, 'TRAC05').squeeze().transpose()
array2d_idepth_iT_modpo4=array2d_idepth_iT_modpo4full[:,array1d_iT1y_iT]
in $\mathrm{ \mu M }$
(ndepths,ndays)=array2d_idepth_iT_modno3.shape
array1d_iT_modnitrate3m=np.zeros(ndays)
x=3
xp=[RC[2],RC[3]]
assert np.all(np.diff(xp)>0), "xp is in decreasing order: %r %r" % (xp[0],xp[1])
for iT in first_year:
fp =[array2d_idepth_iT_modno3[3,iT],array2d_idepth_iT_modno3[2,iT]]
nitrate3m=np.interp(x,xp,fp)
array1d_iT_modnitrate3m[iT]=nitrate3m
(ndepths,ndays)=array2d_idepth_iT_modsioh4.shape
array1d_iT_silicicacid3m=np.zeros(ndays)
x=3
xp=[RC[2],RC[3]]
assert np.all(np.diff(xp)>0), "xp is in decreasing order: %r %r" % (xp[0],xp[1])
for iT in first_year:
fp =[array2d_idepth_iT_modsioh4[3,iT],array2d_idepth_iT_modsioh4[2,iT]]
silica3m=np.interp(x,xp,fp)
array1d_iT_silicicacid3m[iT]=silica3m
in $\mathrm{ \mu M }$
array2d_idepth_iT_modfetfull=netcdf_tools.read_netcdf(carfile, 'TRAC07').squeeze().transpose()
array2d_idepth_iT_modfet=array2d_idepth_iT_modfetfull[:,array1d_iT1y_iT]
indir_exp0='data/DataS8_output_mitgcm/exp0'
chlfile=os.path.join(indir,'chl.0000000000.t001.nc')
array2d_idepth_iT_modchlfull=\
read_mitgcm.get_array2d_idepth_iT_chlfull(chlfile)
array2d_idepth_iT_modchl=array2d_idepth_iT_modchlfull[:,array1d_iT1y_iT]
# get vertically integrated Chl a for one year}
def get_array1d_iT_vintchl(indir,array1d_iT1y_iT,depth_end):
chlfile=os.path.join(indir,'chl.0000000000.t001.nc')
# mg Chl-a m^-3
array2d_idepth_iT_chlfull=\
read_mitgcm.get_array2d_idepth_iT_chlfull(chlfile)
# mg Chl-a m^-2
array1d_iT_vintchlfull\
=vstats.vint(
array2d_idepth_iT_tracer=array2d_idepth_iT_chlfull,
array1d_idepth_delR=drF,
depth_end=depth_end)
array1d_iT_vintchl=array1d_iT_vintchlfull[array1d_iT1y_iT]
return array1d_iT_vintchl
in $\mathrm{ mg\ Chla\ m^{-2} }$
depth_end=40
array1d_iT_modvint40mchl\
=get_array1d_iT_vintchl(indir_exp0,array1d_iT1y_iT,depth_end)
in $\mathrm{ mg\ Chla\ m^{-2} }$
depth_end=100
array1d_iT_modvintchl\
=get_array1d_iT_vintchl(indir_exp0,array1d_iT1y_iT,depth_end)
in $\mathrm{ mg\ Chla\ m^{-2} }$
indir_exp1='data/DataS8_output_mitgcm/exp1'
depth_end=100
array1d_iT_exp1vintchl\
=get_array1d_iT_vintchl(indir_exp1,array1d_iT1y_iT,depth_end)
in $\mathrm{ mg\ Chla\ m^{-2} }$
indir_exp2='data/DataS8_output_mitgcm/exp2'
depth_end=100
array1d_iT_exp2vintchl\
=get_array1d_iT_vintchl(indir_exp2,array1d_iT1y_iT,depth_end)
in $\mathrm{ mg\ Chla\ m^{-2} }$
indir_exp3times0_25='data/DataS8_output_mitgcm/exp3times0_25'
indir_exp3times0_50='data/DataS8_output_mitgcm/exp3times0_50'
indir_exp3times2_00='data/DataS8_output_mitgcm/exp3times2_00'
indir_exp3times4_00='data/DataS8_output_mitgcm/exp3times4_00'
depth_end=100
array1d_iT_exp3times0_25vintchl\
=get_array1d_iT_vintchl(indir_exp3times0_25,array1d_iT1y_iT,depth_end)
array1d_iT_exp3times0_50vintchl\
=get_array1d_iT_vintchl(indir_exp3times0_50,array1d_iT1y_iT,depth_end)
array1d_iT_exp3times2_00vintchl\
=get_array1d_iT_vintchl(indir_exp3times2_00,array1d_iT1y_iT,depth_end)
array1d_iT_exp3times4_00vintchl\
=get_array1d_iT_vintchl(indir_exp3times4_00,array1d_iT1y_iT,depth_end)
tavefile='data/DataS8_output_mitgcm/exp0/tave.0000000000.t001.nc'
array2d_idepth_iT_Ttavefull=netcdf_tools.read_netcdf(tavefile,'Ttave').squeeze().transpose()
array2d_idepth_iT_Ttavefull[-1,:]=np.nan
array2d_idepth_iT_Ttave=array2d_idepth_iT_Ttavefull[:,array1d_iT1y_iT]
minT=np.nanmin(array2d_idepth_iT_Ttave)
maxT=np.nanmax(array2d_idepth_iT_Ttave)
IPython.display.Markdown(
'''The interval of temperature ($^{\circ} C$) is
$[%f,%f].$'''
%(minT,maxT)
)
The interval of temperature ($^{\circ} C$) is $[-1.755564,2.919083].$
The modification of growth rate by temperature (unitless, between 0 and 1) is
$\gamma^T=\tau_T e^{A_T(\frac{1}{T+273.15}-\frac{1}{T_N})}$
where
$\tau_T$ is the normalization factor for temperature function $=0.8$ (unitless)
$A_T$ is a coefficient $=-4000$ K
$T$ is the ambient temperature (K)
$T_N$ is the reference temperature $=20\mathrm{ ^{\circ}C }+273.15\mathrm{ ^{\circ}C }=293.15$ K
tau_T=0.8 # unitless
A_T=-4000 # K
T_N=20+273.15 # K
mingammaT=tau_T*math.exp(A_T*(1/(minT+273.15)-1/T_N))
maxgammaT=tau_T*math.exp(A_T*(1/(maxT+273.15)-1/T_N))
IPython.display.Markdown(
'''The interval of $\gamma^T$ (unitless)
in the simulation is
$[%f,%f].$'''
%(mingammaT,maxgammaT)
)
The interval of $\gamma^T$ (unitless) in the simulation is $[0.267952,0.343909].$
minpcmaxgammaT=pcmax[15-1]*mingammaT
maxpcmaxgammaT=pcmax[15-1]*maxgammaT
IPython.display.Markdown(
'''The interval of $P^C_{max}\ \gamma^T(T)$ for the diatom 6.6 $\mu$m in $d^{-1}$ is
$[%f,%f].$'''
%(minpcmaxgammaT,maxpcmaxgammaT)
)
The interval of $P^C_{max}\ \gamma^T(T)$ for the diatom 6.6 $\mu$m in $d^{-1}$ is $[0.489756,0.628588].$
alphachl=info.f['alphachl'] # mmol C (mg Chl a)^-1 (umol photons m^-2)^-1
alphachl=alphachl[15-1]
For your information, $\alpha^{chl}$ is the same for all numerical phytoplankton types.
We convert the units to compare with literature.
modalphachl_modunits=alphachl # mmol C (mg Chl a)^-1 (umol photons m^-2)^-1
mmolC2molC=1E-3 # mol C (mmol C )^-1
molarmassC=12.011 # g C (mol C)^-1
mgChla2gChla=1E3 # mg Chl a (g Chl a)^-1
s2h=3600 # sh^-1
modalphachl=modalphachl_modunits*mmolC2molC*molarmassC*mgChla2gChla*s2h
IPython.display.Markdown(
r'''$\alpha^{chl} = %.3G
\ g\ C\ (g\ Chl\ a)^{-1}\ h^{-1}\ (\mu mol\ photons\ m^{-2}\ s^{-1})^{-1}$.'''
%(modalphachl)
)
$\alpha^{chl} = 0.0346 \ g\ C\ (g\ Chl\ a)^{-1}\ h^{-1}\ (\mu mol\ photons\ m^{-2}\ s^{-1})^{-1}$.
Maximum photosynthesis rate at saturation irradiance in nutrient replete conditions for numerical diatom 6.6 $\mu$m in $\mathrm{ d^{-1} }$
$P^C_m(j=15) = P^C_{max}(j=15) * \gamma^T(T)$
array2d_idepth_iT_modPCmrep15=pcmax[15-1] \
* tau_T*np.exp(A_T*(1/(array2d_idepth_iT_Ttave+273.15)-1/T_N))
Chlorophyll $a$ for numerical diatom 6.6 $\mu$m ($\mathrm{ mg\ Chl\ a\ m^{-3} }$)
chlfile='data/DataS8_output_mitgcm/exp0/chl.0000000000.t001.nc'
array2d_idepth_iT_diatom7umchlfull\
=netcdf_tools.read_netcdf(chlfile, 'TRAC84').squeeze().transpose()
array2d_idepth_iT_diatom7umchlfull[-1,:]=np.nan
array2d_idepth_iT_diatom7umchl\
=array2d_idepth_iT_diatom7umchlfull[:,array1d_iT1y_iT]
Biomass for numerical diatom 6.6 $\mu$m ($\mathrm{ mmol\ C\ m^{-3} }$)
array2d_idepth_iT_diatom7umbiofull\
=netcdf_tools.read_netcdf(carfile, 'TRAC35').squeeze().transpose()
array2d_idepth_iT_diatom7umbiofull[-1,:]=np.nan
array2d_idepth_iT_diatom7umbio\
=array2d_idepth_iT_diatom7umbiofull[:,array1d_iT1y_iT]
Chl $a$: C ratio for numerical diatom 6.6 $\mu$m ($\mathrm{ mg\ Chl\ {a}\ (mmol\ C)^{-1} }$)
array2d_idepth_iT_modtheta15\
=array2d_idepth_iT_diatom7umchl/array2d_idepth_iT_diatom7umbio
Light saturation parameter for numerical diatom 6.6 $\mu$m ($E_k$, $\mathrm{ \mu mol\ photons\ m^{-2}\ s^{-1} }$)
$E_k= \frac{P^C_m}{\alpha^{chl} \cdot \theta \cdot \frac{86400\ s}{d}}$
where
$E_k$: light saturation parameter ($\mathrm{ \mu mol\ photons\ m^{-2}\ s^{-1} }$).
$P^C_m$: carbon-specific maximum photosynthesis rate ($\mathrm{ d^{-1} }$).
$\alpha^{chl}$: linear initial slope of the Chl a-specific photosynthesis versus irradiance curve ($\mathrm{ mmol\ C\ (mg\ Chl\ a)^{-1}\ (\mu mol\ photons\ m^{-2})^{-1} }$).
$\theta$: Chl a : C ratio ($\mathrm{ mg\ Chl\ a\ (mmol\ C)^{-1} }$).
s2d=86400 # s d^-1
array2d_idepth_iT_modEk15=\
array2d_idepth_iT_modPCmrep15 \
/ (modalphachl_modunits * array2d_idepth_iT_modtheta15 * s2d)
array2d_idepth_iT_modEk15light=array2d_idepth_iT_modEk15.copy()
array2d_idepth_iT_modEk15light[array2d_idepth_iT_meanPARumolEonm2s<5]=np.nan
def make_plots(axs):
box=dict(pad=5,alpha=0)
locs=np.array([0, 31, 59, 90, 120, 151,
181, 212, 243, 273, 304, 334])
labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
# --- ICE CONCENTRATIONS
ax=axs[0]
h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
ax.set_xlim(first_year[0],first_year[-1])
xlims=ax.get_xlim()
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_yticks([0,50,100])
ax.yaxis.set_ticks_position('both')
ax.set_ylabel('Ice conc.\n (%)',bbox=box)
ax.set_ylim(0,100)
ax.grid()
ax.legend(bbox_to_anchor=(0.21,1.00))
plt.text(0.03,0.7,'(a)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# adding the vertical line at 15% of observed ice
ax.axvline(obs_breakup,color='k',linestyle='--')
# --- SNOW AND ICE THICKNESS
ax=axs[1]
h1=ax.fill_between(first_year,array1d_iT_snow_volume,
linestyle='-',color='aliceblue',linewidth=1)
h2=ax.plot(first_year,array1d_iT_snow_volume,
linestyle='-',color='black',lw=1)
h3=ax.plot((snow_thickness_df['doy']-1).to_numpy(),
snow_thickness_df['sample_thickness_m_average'].to_numpy(),
'o',color='black')
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('Thickness\n(m)',bbox=box)
mpl.axis.Axis.set_label_coords(ax.yaxis,-0.061,0.3)
ax.set_ylim(-1.7,1.0)
ax.grid()
ax.text(60,0.2,'snow',
horizontalalignment = 'left',verticalalignment = 'center',
fontsize=20)
h4=ax.fill_between(first_year,-array1d_iT_ice_volume,
linestyle='-',color='dodgerblue',linewidth=1)
h5=ax.plot(first_year,-array1d_iT_ice_volume,
linestyle='-',color='black',lw=1)
h6=ax.plot((ice_thickness_df['doy']-1).to_numpy(),
-ice_thickness_df['sample_thickness_m_average'].to_numpy(),
'o',color='black')
ax.text(60,-0.30,'ice',
horizontalalignment = 'left',verticalalignment = 'center',
fontsize=20)
ax.legend([h2[0],h3[0]],['model','obs'],loc=4,
bbox_to_anchor=(0.90,0.00))
plt.text(0.03,0.85,'(b)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# adding the vertical line at 15% of observed ice
ax.axvline(obs_breakup,linestyle='--',color='black')
# --- sfc PAR
ax=axs[2]
h1=ax.plot(first_year,array1d_iT_modPAR3m,'-',color='black',lw=3)
ax.set_yscale('log')
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('PAR at 3 m\n(mol photons m$^{-2}$ d$^{-1})$',bbox=box)
ax.yaxis.set_ticks_position('both')
ax.set_ylim(1E-2,50)
ax.grid()
plt.text(0.03,0.90,'(c)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
h2=ax.plot(first_year,array1d_iT_obsPAR3m,
'o',color='black')
ax.legend([h1[0],h2[0]],['model','obs'],loc=0, ncol=1)
# adding the vertical line at 15% of observed ice
ax.axvline(obs_breakup,color='k',linestyle='--')
# --- sfc NO3 + sfc SiO2
ax=axs[3]
h1=ax.plot(first_year,array1d_iT_modnitrate3m,'-',color='blue',lw=3)
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('Nutrients at 3 m\n($\mathrm{ \mu M }$)',bbox=box)
ax.yaxis.set_ticks_position('both')
ax.set_ylim(0,10)
ax.grid()
plt.text(0.03,0.90,'(d)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
h2=ax.plot(first_year,array1d_iT_silicicacid3m ,'-',color='red', lw=3)
h3=ax.plot(first_year,array1d_iT_obsnitrate3m,
'o',color='blue')
h4=ax.plot(first_year,array1d_iT_obssilica3m ,
's',color='red')
ax.legend([h1[0],h3[0],
h2[0],h4[0]],
['model $\mathrm{ NO_3 }$','obs $\mathrm{ NO_3 }$',
'model $\mathrm{ Si{(OH)}_4 }$','obs $\mathrm{ Si{(OH)}_4 }$'],
loc=1, ncol=2)
# adding the vertical line at 15% of observed ice
ax.axvline(obs_breakup,color='k',linestyle='--')
# --- CHLOROPHYLL VINT
ax=axs[4]
h1=ax.plot(first_year,array1d_iT_modvint40mchl,
'--',lw=3,color='black',label='model (0-40m)')
h2=ax.plot(first_year,array1d_iT_modvintchl,
'-',color='black',lw=3,label='model (0-100m)')
h3=ax.plot(first_year,array1d_iT_obsvint40mchlHPLC,
'+',color='black',label='obs (0-40m)')
h4=ax.plot(first_year,array1d_iT_obsvintchlHPLC,
'o',color='black',label='obs (0-100m)')
ax.set_yscale('log')
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('Vertically integrated Chl $a$\n($\mathrm{ mg\ Chl\ m^{-2} }$)',
bbox=box)
ax.yaxis.set_ticks_position('both')
ax.set_ylim(0.2,300)
ax.grid()
plt.text(0.03,0.90,'(e)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
ax.legend(bbox_to_anchor=(0.31,1.00))
# adding the vertical line at 15% of observed ice
ax.axvline(obs_breakup,color='k',linestyle='--')
# -- ISOLUME
ax=axs[5]
h1=ax.plot(first_year,
array1d_iT_modisolume,
'-',
lw=3,
color='black',
label='model')
h2=ax.plot(first_year,
array1d_iT_obsisolume,
'o',
color='black',
label='obs')
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('Depth of $\mathrm{ z_{0.415}\ isolume }$\n(m)',
bbox=box)
ax.yaxis.set_ticks_position('both')
ax.set_ylim(0,80)
ax.invert_yaxis()
ax.legend(loc=3)
ax.grid()
plt.text(0.03,0.90,'(f)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# adding the vertical line at 15% of ice
ax.axvline(obs_breakup,color='k',linestyle='--')
# --- ADDITIONAL X-AXIS
ax=axs[6]
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels(labels)
ax.axes.get_yaxis().set_visible(False)
ax.spines['top'].set_visible(False)
# --- POSITION
axs[0].set_position( [0.13,0.95 ,0.85,0.03])
axs[1].set_position( [0.13,0.88 ,0.85,0.06])
axs[2].set_position( [0.13,0.67 ,0.85,0.20])
axs[3].set_position( [0.13,0.46 ,0.85,0.20])
axs[4].set_position( [0.13,0.25 ,0.85,0.20])
axs[5].set_position( [0.13,0.04 ,0.85,0.20])
axs[6].set_position( [0.13,0.02 ,0.85,0.01])
with plt.style.context('mplstyles/obs_and_model.mplstyle'):
# Plot
fig,axs=plt.subplots(7,1,sharex=False,figsize=(16, 24))
make_plots(axs)
fig.align_ylabels(axs[:])
# --- SAVE
plt.savefig('figures_progress/qik.png',dpi=1000)
type_id=np.arange(1,52)
pos_Prochlorococcus=np.arange( 0, 1)
pos_Synechococcus =np.arange( 1, 2)
pos_smalleuk =np.arange( 2, 4)
pos_cocco =np.arange( 4, 9)
pos_diazo =np.arange( 9,13)
pos_Trichodesmium =np.arange(13,14)
pos_diatoms =np.arange(14,23)
pos_mixodino =np.arange(23,33)
pos_zoo =np.arange(33,49)
def make_plots(axs):
box=dict(pad=5,alpha=0)
# --- EQUIVALENT SPHERICAL DIAMETER (ESD)
ax=axs[0]
h3=ax.plot(type_id[pos_smalleuk]-2, esd[pos_smalleuk],'gh',markersize=12)
h4=ax.plot(type_id[pos_cocco]-2, esd[pos_cocco], 'b+')
h7=ax.plot(type_id[pos_diatoms]-2-5, esd[pos_diatoms], 'ro')
h8=ax.plot(type_id[pos_mixodino]-2-5,esd[pos_mixodino],'md')
h9=ax.plot(type_id[pos_zoo]-2-5, esd[pos_zoo], 'k*')
ax.set_xticklabels([])
ax.xaxis.set_major_locator(mticker.MultipleLocator(10))
ax.xaxis.set_minor_locator(mticker.MultipleLocator(1))
ax.set_xlim([0,43])
ax.set_ylabel('Equivalent spherical diameter ($\mathrm{ \mu m }$)')
ax.set_yscale('log')
ax.grid()
plt.text(0.03,0.92,'(a)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# --- MAXIMUM GROWTH RATE (P^C_MAX)
ax=axs[1]
h3=ax.plot(type_id[pos_smalleuk]-2, pcmax[pos_smalleuk],'gh',markersize=12)
h4=ax.plot(type_id[pos_cocco]-2, pcmax[pos_cocco], 'b+')
h7=ax.plot(type_id[pos_diatoms]-2-5, pcmax[pos_diatoms], 'ro')
h8=ax.plot(type_id[pos_mixodino]-2-5,pcmax[pos_mixodino],'md')
h9=ax.plot(type_id[pos_zoo]-2-5, pcmax[pos_zoo], 'k*')
ax.set_xticklabels([])
ax.xaxis.set_major_locator(mticker.MultipleLocator(10))
ax.xaxis.set_minor_locator(mticker.MultipleLocator(1))
ax.set_xlim([0,43])
ax.set_ylabel('Maximum growth rate ($\mathrm{ d^{-1} }$)')
ax.set_yscale('log')
ax.yaxis.set_major_locator(plt.FixedLocator([0.5,1,1.5,2]))
ax.yaxis.set_minor_locator(plt.FixedLocator([0.6,0.7,0.8,0.9,
1.1,1.2,1.3,1.4,1.6,1.7,1.8,1.9]))
ax.yaxis.set_major_formatter(plt.FixedFormatter([0.5,1,1.5,2]))
ax.yaxis.set_minor_formatter(plt.FixedFormatter(['']))
ax.set_ylim([0.5,2])
ax.grid()
plt.text(0.03,0.92,'(b)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# --- HALF SATURATION FOR GROWTH ON NITRATE (k_NO_3)
ax=axs[2]
h3=ax.plot(type_id[pos_smalleuk]-2, ksatno3[pos_smalleuk],'gh',markersize=12)
h4=ax.plot(type_id[pos_cocco]-2, ksatno3[pos_cocco], 'b+')
h7=ax.plot(type_id[pos_diatoms]-2-5, ksatno3[pos_diatoms], 'ro')
h8=ax.plot(type_id[pos_mixodino]-2-5,ksatno3[pos_mixodino],'md')
h9=ax.plot(type_id[pos_zoo]-2-5, ksatno3[pos_zoo], 'k*')
ax.set_xticklabels([])
ax.xaxis.set_major_locator(mticker.MultipleLocator(10))
ax.xaxis.set_minor_locator(mticker.MultipleLocator(1))
ax.set_xlim([0,43])
ax.set_ylabel('$\mathrm{ k_{NO_3} }$ ($\mathrm{ mmol\ N\ m^{-3} }$)')
ax.set_yscale('log')
ax.set_ylim([0.01,10])
ax.grid()
plt.text(0.03,0.92,'(c)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# --- MAXIMUM GRAZING RATE (G_MAX)
ax=axs[3]
h3=ax.plot(type_id[pos_smalleuk]-2, gmax[pos_smalleuk],'gh',markersize=12,
label='picophytoplankton')
h4=ax.plot(type_id[pos_cocco]-2, gmax[pos_cocco], 'b+',
label='other nanophytoplankton')
h7=ax.plot(type_id[pos_diatoms]-2-5, gmax[pos_diatoms], 'ro',
label='diatoms')
h8=ax.plot(type_id[pos_mixodino]-2-5,gmax[pos_mixodino],'md',
label='mixotrophic dinoflagellates')
h9=ax.plot(type_id[pos_zoo]-2-5, gmax[pos_zoo], 'k*',
label='microzooplankton')
ax.set_xlabel('Planktonic type ID',fontsize=23)
ax.xaxis.set_major_locator(mticker.MultipleLocator(10))
ax.xaxis.set_minor_locator(mticker.MultipleLocator(1))
ax.set_xlim([0,43])
ax.set_ylabel('Maximum grazing rate ($\mathrm{ d^{-1} }$)')
ax.set_yscale('log')
ax.yaxis.set_major_locator(plt.FixedLocator([0.5,1,5,10]))
ax.yaxis.set_minor_locator(plt.FixedLocator([0.6,0.7,0.8,0.9,2,3,4,6,7,8,9,20]))
ax.yaxis.set_major_formatter(plt.FixedFormatter([0.5,1,5,10]))
ax.set_ylim([0.5,20])
ax.grid()
ax.legend(loc=3)
plt.text(0.03,0.92,'(d)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# --- POSITION
axs[0].set_position( [0.10,0.76 ,0.85,0.21])
axs[1].set_position( [0.10,0.53 ,0.85,0.21])
axs[2].set_position( [0.10,0.30 ,0.85,0.21])
axs[3].set_position( [0.10,0.07 ,0.85,0.21])
with plt.style.context('mplstyles/gud_traits.mplstyle'):
# Plot
fig,axs=plt.subplots(4,1,figsize=(16, 16))
make_plots(axs)
fig.align_ylabels(axs[:])
# --- SAVE
plt.savefig('figures_progress/traits.png',dpi=1000)
def make_plots(axs):
locs=np.array([0, 31, 59, 90, 120, 151,
181, 212, 243, 273, 304, 334])
labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
# --- NITRATE
ax=axs[0]
h=ax.pcolormesh(first_year366,
RF_above151,
array2d_idepth_iT_modno3[0:(RF_above151.size)-1,:],
cmap='viridis',
vmin=0,
vmax=6)
ax.scatter(nutrients_df['doy']-1,
nutrients_df['depth'],
c=nutrients_df['no3'],
cmap='viridis',
vmin=0,
vmax=6,
s=100,
edgecolor='k')
ax.set_xticklabels([])
ax.set_ylabel('Depth (m)')
ax.set_xlim(first_year366[0],first_year[-1])
ax.set_ylim(150,0)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
cbaxes = fig.add_axes([0.88, 0.83, 0.01, 0.14])
cbar=plt.colorbar(h,
cax = cbaxes)
cbar.set_label('(mmol N $\mathrm{ m^{-3} }$)')
plt.text(0.50,1.10,'(a) nitrate',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes)
# --- SILICIC ACID
ax=axs[1]
h=ax.pcolormesh(first_year366,
RF_above151,
array2d_idepth_iT_modsioh4[0:(RF_above151.size)-1,:],
cmap='viridis',
vmin=0,
vmax=10)
ax.scatter(nutrients_df['doy']-1,
nutrients_df['depth'],
c=nutrients_df['sioh4'],
cmap='viridis',
vmin=0,
vmax=10,
s=100,
edgecolor='k')
ax.set_xticklabels([])
ax.set_ylabel('Depth (m)')
ax.set_xlim(first_year366[0],first_year[-1])
ax.set_ylim(150,0)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
cbaxes = fig.add_axes([0.88, 0.63, 0.01, 0.14])
cbar=plt.colorbar(h,
cax = cbaxes)
cbar.set_label('(mmol Si$\mathrm{ {(OH)}_4 }$ $\mathrm{ m^{-3} }$)')
plt.text(0.50,1.10,'(b) silicic acid',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes)
# --- PHOSPHATE
ax=axs[2]
h=ax.pcolormesh(first_year366,
RF_above151,
array2d_idepth_iT_modpo4[0:(RF_above151.size)-1,:],
cmap='viridis',
vmin=0,
vmax=1)
ax.scatter(nutrients_df['doy']-1,
nutrients_df['depth'],
c=nutrients_df['po4'],
cmap='viridis',
vmin=0,
vmax=1,
s=100,
edgecolor='k')
ax.set_xticklabels([])
ax.set_ylabel('Depth (m)')
ax.set_xlim(first_year366[0],first_year[-1])
ax.set_ylim(150,0)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
cbaxes = fig.add_axes([0.88, 0.43, 0.01, 0.14])
cbar=plt.colorbar(h,
cax = cbaxes)
cbar.set_label('(mmol $\mathrm{ PO_4 }$ $\mathrm{ m^{-3} }$)')
plt.text(0.50,1.10,'(c) phosphate',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes)
# --- IRON
ax=axs[3]
h=ax.pcolormesh(first_year366,
RF_above151,
array2d_idepth_iT_modfet[0:(RF_above151.size)-1,:],
cmap='viridis',
vmin=0,
vmax=1.5E-3)
ax.scatter(fet_df['doy']-1,
fet_df['depth'],
c=fet_df['Fe_D_CONC_BOTTLE_mmol_m3'],
cmap='viridis',
vmin=0,
vmax=1.5E-3,
s=100,
edgecolor='k')
ax.set_xticklabels([])
ax.set_ylabel('Depth (m)')
ax.set_xlim(first_year366[0],first_year[-1])
ax.set_ylim(150,0)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
cbaxes = fig.add_axes([0.88, 0.23, 0.01, 0.14])
cbar=plt.colorbar(h,
cax = cbaxes)
cbar.set_label('(mmol Fe $\mathrm{ m^{-3} }$)')
plt.text(0.50,1.10,'(d) total dissolved inorganic iron',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes)
# --- CHLOROPHYLL A
ax=axs[4]
h=ax.pcolormesh(first_year366,
RF_above151,
array2d_idepth_iT_modchl[0:(RF_above151.size)-1,:],
cmap='viridis',
norm=mpl.colors.LogNorm(vmin=1E-2, vmax=10))
ax.scatter(chlHPLC_df['doy']-1,
chlHPLC_df['depth'],
c=chlHPLC_df['chlHPLC'],
cmap='viridis',
norm=mpl.colors.LogNorm(vmin=1E-2, vmax=10),
s=100,
edgecolor='k')
ax.set_xticklabels([])
ax.set_ylabel('Depth (m)')
ax.set_xlim(first_year366[0],first_year[-1])
ax.set_ylim(150,0)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
cbaxes = fig.add_axes([0.88, 0.03, 0.01, 0.14])
cbar=plt.colorbar(h,
cax = cbaxes)
cbar.set_label('(mg Chl $a$ $\mathrm{ m^{-3} }$)')
plt.text(0.50,1.10,'(e) chlorophyll $a$',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes)
# --- POSITION
axs[0].set_position( [0.07,0.83 ,0.79,0.14])
axs[1].set_position( [0.07,0.63 ,0.79,0.14])
axs[2].set_position( [0.07,0.43 ,0.79,0.14])
axs[3].set_position( [0.07,0.23 ,0.79,0.14])
axs[4].set_position( [0.07,0.03 ,0.79,0.14])
with plt.style.context('mplstyles/heatmaps.mplstyle'):
# Plot
fig,axs=plt.subplots(5,1,sharex=False,figsize=(16, 24))
make_plots(axs)
# --- SAVE
plt.savefig('figures_progress/conc.png',dpi=1000)
obsPCmax_df=obsPCmax_df.copy() # to avoid a SettingWithCopyWarning
obsPCmax_df.sort_values(by=['species'],inplace=True)
obsalphachl_df=obsalphachl_df.copy() # to avoid a SettingWithCopyWarning
obsalphachl_df.sort_values(by=['species'],inplace=True)
obsEk_df=obsEk_df.copy() # to avoid a SettingWithCopyWarning
obsEk_df.sort_values(by=['species'],inplace=True)
array1d_i_modEk15light=array2d_idepth_iT_modEk15light.flatten()
pos=~np.isnan(array1d_i_modEk15light)
array1d_i_modEk15light=array1d_i_modEk15light[pos]
def make_plots(axs):
# --- MAXIMUM PHOTOSYNTHESIS RATE AT SATURATION IRRADIANCE
# --- IN NUTRIENT REPLETE CONDITIONS
# --- (P^C_max * gamma^T)
ax=axs[0,0]
xlim=(0,1.5)
h1=ax.scatter(obsPCmax_df.value,obsPCmax_df.species,
color='b',label='obs')
h2=ax.axhline('',
xmin=minpcmaxgammaT/xlim[1],
xmax=maxpcmaxgammaT/xlim[1],
color='red',lw=10)
h3=ax.axhline('foo',
lw=0)
# to have the ytick labels in alphabetical order from top to bottom
ax.invert_yaxis()
# in italics
# FULL GENUS
ax.set_yticklabels(['Chaetoceros furcellatus',
'Chaetoceros neogracile',
'Fragilariopsis cylindrus',
'Thalassiosira gravida 1',
'Thalassiosira gravida 2',
'Thalassiosira nordenskioeldii',
'Numerical diatom 6.6 μm'
],
fontstyle='italic',
fontweight='bold')
ax.set_xlim(xlim)
ax.set_xlabel('$P^C_{max} \gamma^T$\n($\mathrm{ d^{-1} }$)')
ax.grid()
plt.text(0.10,0.88,'(a)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
ax.legend([h1,h2],['obs','model'],ncol=2,loc='lower center')
# --- NOTHING
axs[1,0].axis('off')
# --- LINEAR INITIAL SLOPE OF THE CHL A-SPECIFIC PHOTOSYNTHESIS
# --- VERSUS IRRADIANCE CURVE
# --- (ALPHA^CHL)
ax=axs[0,1]
xlim=(0,0.05)
h1=ax.scatter(obsalphachl_df.value,
obsalphachl_df.species,
c='b',label='obs')
h2=ax.scatter(modalphachl,
'',
color='red')
h3=ax.axhline('foo',
lw=0)
# to have the ytick labels in alphabetical order from top to bottom
ax.invert_yaxis()
ax.set_yticklabels([])
ax.set_xlim(xlim)
ax.set_xlabel('''$\\alpha^{chl}$\n\
$\mathrm{ (g\ C\ (g\ Chl)^{-1}\ h^{-1} }$\n\
$\mathrm{ (\mu mol\ photons\ m^{-2}\ s^{-1})^{-1}) }$''')
ax.grid()
plt.text(0.10,0.88,'(b)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# --- NOTHING
axs[1,1].axis('off')
# --- LIGHT SATURATION PARAMETER
# --- (E_k)
ax=axs[0,2]
xlim=(0,250)
h1=ax.scatter(obsEk_df.value,
obsEk_df.species,
c='b',label='obs')
h3=ax.axhline('foo',
lw=0)
h4=ax.axhline('foo2',
lw=0)
# to have the ytick labels in alphabetical order from top to bottom
ax.invert_yaxis()
ax.set_yticklabels([])
ax.set_xlabel('$E_k$\n$\mathrm{ (\mu mol\ photons\ m^{-2}\ s^{-1}) }$')
ax.set_xlim(xlim)
ax.grid()
plt.text(0.10,0.88,'(c)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
ax=axs[1,2]
boxplot=ax.boxplot(array1d_i_modEk15light,
vert=False,
widths=0.75,
boxprops=dict(color='red'),
flierprops=dict(marker='o',
markerfacecolor='red',
markeredgecolor='none',
markersize=1),
medianprops=dict(color='red'),
capprops=dict(color='red'),
whiskerprops=dict(color='red'))
ax.set_xlim(xlim)
ax.set_yticklabels([])
ax.set_xticklabels([])
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['bottom'].set_visible(False)
# ax.spines['left'].set_visible(False)
ax.tick_params(length=0)
ax.grid()
# --- POSITION
# FULL GENUS
axs[0,0].set_position( [0.24,0.21 ,0.22,0.77])
axs[0,1].set_position( [0.50,0.21 ,0.22,0.77])
axs[0,2].set_position( [0.76,0.21 ,0.22,0.77])
axs[1,2].set_position( [0.76,0.31 ,0.22,0.07])
return boxplot
with plt.style.context('mplstyles/lacour17.mplstyle'):
# Plot
fig,axs=plt.subplots(2,3,figsize=(16, 8))
boxplot=make_plots(axs)
# --- SAVE
plt.savefig('figures_progress/lacour17.png',dpi=1000)
Where IQR is the interquartile range (Q3-Q1), the upper whisker will extend to last datum less than Q3+whis*IQR (https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.boxplot.html).
upper_whisker=[item.get_xdata()[1] for item in boxplot['whiskers']][1]
fliers=array1d_i_modEk15light[array1d_i_modEk15light>upper_whisker]
nb_fliers=fliers.size
nb_total=array1d_i_modEk15light.size
Number of outliers
nb_fliers
311
Number of data points
nb_total
3409
data={'exp0x':[0, 163, 198,228,300,328,364],
'exp0y':[1, 1, 100,100, 10, 1, 1],
'exp1x':[0, 201, 225,255,300,328,364],
'exp1y':[1, 1, 100,100, 10, 1, 1],
'exp2x':[0, 163, 228,258,300,328,364],
'exp2y':[0.001,0.001,100,100, 10, 1, 1],
'exp3x':[0, 163, 181,210,300,325,364],
'exp3y':[1, 1, 10, 10, 3, 1, 1]}
df=pd.DataFrame(data)
# df
def make_plots(axs):
locs=np.array([0, 31, 59, 90, 120, 151,
181, 212, 243, 273, 304, 334])
labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
xlims=(first_year[0],first_year[-1])
xmin=xlims[0]
xmax=xlims[1]
# --- EXP-0: ICE CONCENTRATION
ax=axs[0,0]
h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
ax.set_xlim(xlims)
ymin=0
ymax=100
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_yticks([0,50,100])
ax.yaxis.set_ticks_position('both')
ax.set_ylabel('Ice conc.\n (%)')
ax.set_ylim(0,100)
ax.grid()
ax.legend(loc='lower left',bbox_to_anchor=(0.08,-0.15))
ax.axvline(obs_breakup,color='k',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
plt.text(0.03,0.60,'(a)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
plt.text(0.5,1.5,'EXP-0: default',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- EXP-1: ICE CONCENTRATION
ax=axs[0,1]
h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
ax.set_xlim(xlims)
ymin=0
ymax=100
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_yticks([0,50,100])
ax.yaxis.set_ticks_position('both')
ax.set_ylabel('Ice conc.\n (%)')
ax.set_ylim(ymin,ymax)
ax.grid()
ax.legend(loc='lower left',bbox_to_anchor=(0.08,-0.15))
ax.axvline(obs_breakup,color='k',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
plt.text(0.03,0.60,'(b)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
plt.text(0.5,1.5,'EXP-1: no light under sea ice',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- EXP-0: CHLOROPHYLL VINT
ax=axs[1,0]
h1=ax.plot(first_year,array1d_iT_modvintchl,
'-',color='black',lw=3,label='model (0-100m)')
h2=ax.plot(first_year,array1d_iT_obsvintchlHPLC,
'o',color='black',label='obs (0-100m)')
ax.set_yscale('log')
ymin=0.2
ymax=300
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('Vertically integrated Chl $a$\n($\mathrm{ mg\ Chl\ m^{-2} }$)')
ax.yaxis.set_ticks_position('both')
ax.set_ylim(ymin,ymax)
ax.grid()
ax.legend(loc='lower left',bbox_to_anchor=(0.08,0.75))
ax.axvline(obs_breakup,color='k',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
plt.text(0.03,0.90,'(c)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- EXP-1: CHLOROPHYLL VINT
ax=axs[1,1]
h1=ax.plot(first_year,array1d_iT_exp1vintchl,
'-.',color='red',lw=3,label='model (0-100m)')
h2=ax.plot(first_year,array1d_iT_obsvintchlHPLC,
'o',color='black',label='obs (0-100m)')
ax.set_yscale('log')
ymin=0.2
ymax=300
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('Vertically integrated Chl $a$\n($\mathrm{ mg\ Chl\ m^{-2} }$)')
ax.yaxis.set_ticks_position('both')
ax.set_ylim(ymin,ymax)
ax.grid()
ax.legend(loc='lower left',bbox_to_anchor=(0.08,0.75))
ax.axvline(obs_breakup,color='k',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
plt.text(0.03,0.90,'(d)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- EXP-0: ADDITIONAL X-AXIS
ax=axs[2,0]
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels(labels)
ax.axes.get_yaxis().set_visible(False)
ax.spines['top'].set_visible(False)
# --- EXP-1: ADDITIONAL X-AXIS
ax=axs[2,1]
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels(labels)
ax.axes.get_yaxis().set_visible(False)
ax.spines['top'].set_visible(False)
# --- EXP-0: SCHEMATIC
ax=axs[3,0]
h0=ax.plot(np.array(df.exp0x),np.array(df.exp0y),
'-',color='black',lw=3,label='EXP-0')
ax.set_xlabel('Time')
ax.xaxis.set_label_coords(0.05,-0.05)
ax.set_xlim(xlims)
xticks=([df.exp0x.iloc[1],df.exp0x.iloc[2]])
xticklabels=[r'$\leftarrow t_I \rightarrow$',
r'$\leftarrow t_P \rightarrow$']
ax.set_xticks(xticks)
ax.set_xticklabels(xticklabels)
ax.set_yscale('log')
ax.set_ylabel('Log10(Chlorophyll $a$)')
ax.yaxis.set_label_coords(-0.075,0.50)
ymin=1E-4
ymax=300
ax.set_ylim(ymin,ymax)
ax.tick_params(axis='y',which='both',left=False,labelleft=False)
ax.plot(df.exp0x.iloc[1],df.exp0y.iloc[1],
'o',color='black')
ax.text(df.exp0x.iloc[1]-6,df.exp0y.iloc[1]+1,'I')
ax.plot(df.exp0x.iloc[2],df.exp0y.iloc[2],
'o',color='black')
ax.text(df.exp0x.iloc[2]-6,df.exp0y.iloc[2]+20,'P')
ax.axvline(df.exp0x.iloc[1],
ymax=( math.log10(df.exp0y.iloc[1])-math.log10(ymin) ) \
/ ( math.log10(ymax) -math.log10(ymin) ),
color='black',linestyle='--')
ax.axvline(df.exp0x.iloc[2],
ymax=( math.log10(df.exp0y.iloc[2])-math.log10(ymin) ) \
/ ( math.log10(ymax) -math.log10(ymin) ),
color='black',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
ax.legend(loc='upper right')
plt.text(0.03,0.90,'(e)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- EXP-1: SCHEMATIC
ax=axs[3,1]
h0=ax.plot(np.array(df.exp1x),np.array(df.exp1y),
'-.',color='red',lw=3,label='EXP-1')
ax.set_xlabel('Time')
ax.xaxis.set_label_coords(0.05,-0.05)
ax.set_xlim(xlims)
xticks=([df.exp1x.iloc[1],df.exp1x.iloc[2]])
xticklabels=[r'$\leftarrow t_I \rightarrow$',
r'$\leftarrow t_P \rightarrow$']
ax.set_xticks(xticks)
ax.set_xticklabels(xticklabels)
ax.set_yscale('log')
ax.set_ylabel('Log10(Chlorophyll $a$)')
ax.yaxis.set_label_coords(-0.075,0.50)
ymin=1E-4
ymax=300
ax.set_ylim(ymin,ymax)
ax.tick_params(axis='y',which='both',left=False,labelleft=False)
ax.plot(df.exp1x.iloc[1],df.exp1y.iloc[1],
'o',color='red')
ax.text(df.exp1x.iloc[1]-6,df.exp1y.iloc[1]+1,'I\'',color='red')
ax.plot(df.exp1x.iloc[2],df.exp1y.iloc[2],
'o',color='red')
ax.text(df.exp1x.iloc[2]-6,df.exp1y.iloc[2]+20,'P\'',color='red')
ax.axvline(df.exp1x.iloc[1],
ymax=( math.log10(df.exp1y.iloc[1])-math.log10(ymin) ) \
/ ( math.log10(ymax) -math.log10(ymin) ),
color='black',linestyle='--')
ax.axvline(df.exp1x.iloc[2],
ymax=( math.log10(df.exp1y.iloc[2])-math.log10(ymin) ) \
/ ( math.log10(ymax) -math.log10(ymin) ),
color='black',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
ax.legend(loc='upper right')
plt.text(0.03,0.90,'(f)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- POSITION
axs[0,0].set_position( [0.05,0.90 ,0.43,0.05])
axs[0,1].set_position( [0.54,0.90 ,0.43,0.05])
axs[1,0].set_position( [0.05,0.49 ,0.43,0.38])
axs[1,1].set_position( [0.54,0.49 ,0.43,0.38])
axs[2,0].set_position( [0.05,0.47 ,0.43,0.01])
axs[2,1].set_position( [0.54,0.47 ,0.43,0.01])
axs[3,0].set_position( [0.05,0.04 ,0.43,0.38])
axs[3,1].set_position( [0.54,0.04 ,0.43,0.38])
with plt.style.context('mplstyles/exp.mplstyle'):
# Plot
fig,axs=plt.subplots(4,2,figsize=(31.6,12.3))
make_plots(axs)
# --- SAVE
plt.savefig('figures_progress/exp1.png',dpi=1000)
def make_plots(axs):
locs=np.array([0, 31, 59, 90, 120, 151,
181, 212, 243, 273, 304, 334])
labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
xlims=(first_year[0],first_year[-1])
xmin=xlims[0]
xmax=xlims[1]
# --- EXP-0: ICE CONCENTRATION
ax=axs[0,0]
h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
ax.set_xlim(xlims)
ymin=0
ymax=100
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_yticks([0,50,100])
ax.yaxis.set_ticks_position('both')
ax.set_ylabel('Ice conc.\n (%)')
ax.set_ylim(0,100)
ax.grid()
ax.legend(loc='lower left',bbox_to_anchor=(0.08,-0.15))
ax.axvline(obs_breakup,color='k',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
plt.text(0.03,0.60,'(a)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
plt.text(0.5,1.5,'EXP-0: default',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- EXP-2: ICE CONCENTRATION
ax=axs[0,1]
h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
ax.set_xlim(xlims)
ymin=0
ymax=100
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_yticks([0,50,100])
ax.yaxis.set_ticks_position('both')
ax.set_ylabel('Ice conc.\n (%)')
ax.set_ylim(ymin,ymax)
ax.grid()
ax.legend(loc='lower left',bbox_to_anchor=(0.08,-0.15))
ax.axvline(obs_breakup,color='k',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
plt.text(0.03,0.60,'(b)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
plt.text(0.5,1.5,'EXP-2: no minimum biomass',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- EXP-0: CHLOROPHYLL VINT
ax=axs[1,0]
h1=ax.plot(first_year,array1d_iT_modvintchl,
'-',color='black',lw=3,label='model (0-100m)')
h2=ax.plot(first_year,array1d_iT_obsvintchlHPLC,
'o',color='black',label='obs (0-100m)')
ax.set_yscale('log')
ymin=1E-4
ymax=300
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('Vertically integrated Chl $a$\n($\mathrm{ mg\ Chl\ m^{-2} }$)')
ax.yaxis.set_ticks_position('both')
ax.set_ylim(ymin,ymax)
ax.grid()
ax.legend(loc='lower left',bbox_to_anchor=(0.08,0.75))
ax.axvline(obs_breakup,color='k',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
plt.text(0.03,0.90,'(c)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- EXP-2: CHLOROPHYLL VINT
ax=axs[1,1]
h1=ax.plot(first_year,array1d_iT_exp2vintchl,
'--',color='blue',lw=3,label='model (0-100m)')
h2=ax.plot(first_year,array1d_iT_obsvintchlHPLC,
'o',color='black',label='obs (0-100m)')
ax.set_yscale('log')
ymin=1E-4
ymax=300
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('Vertically integrated Chl $a$\n($\mathrm{ mg\ Chl\ m^{-2} }$)')
ax.yaxis.set_ticks_position('both')
ax.set_ylim(ymin,ymax)
ax.grid()
ax.legend(loc='lower left',bbox_to_anchor=(0.08,0.75))
ax.axvline(obs_breakup,color='k',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
plt.text(0.03,0.90,'(d)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- EXP-0: ADDITIONAL X-AXIS
ax=axs[2,0]
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels(labels)
ax.axes.get_yaxis().set_visible(False)
ax.spines['top'].set_visible(False)
# --- EXP-2: ADDITIONAL X-AXIS
ax=axs[2,1]
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels(labels)
ax.axes.get_yaxis().set_visible(False)
ax.spines['top'].set_visible(False)
# --- EXP-0: SCHEMATIC
ax=axs[3,0]
h0=ax.plot(np.array(df.exp0x),np.array(df.exp0y),
'-',color='black',lw=3,label='EXP-0')
ax.set_xlabel('Time')
ax.xaxis.set_label_coords(0.05,-0.05)
ax.set_xlim(xlims)
xticks=([df.exp0x.iloc[2]])
xticklabels=[r'$\leftarrow t_P \rightarrow$']
ax.set_xticks(xticks)
ax.set_xticklabels(xticklabels)
ax.set_yscale('log')
ax.set_ylabel('Log10(Chlorophyll $a$)')
ax.yaxis.set_label_coords(-0.075,0.50)
ymin=1E-4
ymax=300
ax.set_ylim(ymin,ymax)
ax.tick_params(axis='y',which='minor',left=False,labelleft=False)
yticks=([df.exp0y.iloc[0]])
yticklabels=[r'$\uparrow$'
'\n'
r'$y_I$'
'\n'
r'$\downarrow$']
ax.set_yticks(yticks)
ax.set_yticklabels(yticklabels)
ax.plot(df.exp0x.iloc[1],df.exp0y.iloc[1],
'o',color='black')
ax.text(df.exp0x.iloc[1]-6,df.exp0y.iloc[1]+1,'I')
ax.plot(df.exp0x.iloc[2],df.exp0y.iloc[2],
'o',color='black')
ax.text(df.exp0x.iloc[2]-6,df.exp0y.iloc[2]+20,'P')
ax.axvline(df.exp0x.iloc[2],
ymax=( math.log10(df.exp0y.iloc[2])-math.log10(ymin) ) \
/ ( math.log10(ymax) -math.log10(ymin) ),
color='black',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
ax.legend(loc='upper right')
plt.text(0.03,0.90,'(e)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- EXP-2: SCHEMATIC
ax=axs[3,1]
h0=ax.plot(np.array(df.exp2x),np.array(df.exp2y),
'--',color='blue',lw=3,label='EXP-1')
ax.set_xlabel('Time')
ax.xaxis.set_label_coords(0.05,-0.05)
ax.set_xlim(xlims)
xticks=([df.exp2x.iloc[2]])
xticklabels=[r'$\leftarrow t_P \rightarrow$']
ax.set_xticks(xticks)
ax.set_xticklabels(xticklabels)
ax.set_yscale('log')
ax.set_ylabel('Log10(Chlorophyll $a$)')
ax.yaxis.set_label_coords(-0.075,0.50)
ymin=1E-4
ymax=300
ax.set_ylim(ymin,ymax)
ax.tick_params(axis='y',which='minor',left=False,labelleft=False)
yticks=([df.exp2y.iloc[0]])
yticklabels=[r'$\uparrow$'
'\n'
r'$y_I$'
'\n'
r'$\downarrow$']
ax.set_yticks(yticks)
ax.set_yticklabels(yticklabels)
ax.plot(df.exp2x.iloc[1],df.exp2y.iloc[1],
'o',color='blue')
ax.text(df.exp2x.iloc[1]-6,df.exp2y.iloc[1]+1E-3,'I\'\'',color='blue')
ax.plot(df.exp2x.iloc[2],df.exp2y.iloc[2],
'o',color='blue')
ax.text(df.exp2x.iloc[2]-6,df.exp2y.iloc[2]+20,'P\'\'',color='blue')
ax.axvline(df.exp2x.iloc[2],
ymax=( math.log10(df.exp2y.iloc[2])-math.log10(ymin) ) \
/ ( math.log10(ymax) -math.log10(ymin) ),
color='black',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
ax.legend(loc='upper right')
plt.text(0.03,0.90,'(f)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- POSITION
axs[0,0].set_position( [0.05,0.90 ,0.43,0.05])
axs[0,1].set_position( [0.54,0.90 ,0.43,0.05])
axs[1,0].set_position( [0.05,0.49 ,0.43,0.38])
axs[1,1].set_position( [0.54,0.49 ,0.43,0.38])
axs[2,0].set_position( [0.05,0.47 ,0.43,0.01])
axs[2,1].set_position( [0.54,0.47 ,0.43,0.01])
axs[3,0].set_position( [0.05,0.04 ,0.43,0.38])
axs[3,1].set_position( [0.54,0.04 ,0.43,0.38])
with plt.style.context('mplstyles/exp.mplstyle'):
# Plot
fig,axs=plt.subplots(4,2,figsize=(31.6,12.3))
make_plots(axs)
# --- SAVE
plt.savefig('figures_progress/exp2.png',dpi=1000)
def make_plots(axs):
locs=np.array([0, 31, 59, 90, 120, 151,
181, 212, 243, 273, 304, 334])
labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
xlims=(first_year[0],first_year[-1])
xmin=xlims[0]
xmax=xlims[1]
# --- EXP-0: ICE CONCENTRATION
ax=axs[0]
h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
ax.set_xlim(xlims)
ymin=0
ymax=100
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_yticks([0,50,100])
ax.yaxis.set_ticks_position('both')
ax.set_ylabel('Ice conc.\n (%)')
ax.set_ylim(0,100)
ax.grid()
ax.legend(loc='lower left',bbox_to_anchor=(0.08,-0.15))
ax.axvline(obs_breakup,color='k',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
plt.text(0.03,0.60,'(a)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
plt.text(0.5,1.5,'EXP-3: differing wintertime nitrate',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- EXP-0: CHLOROPHYLL VINT
ax=axs[1]
h0_25=ax.plot(first_year,array1d_iT_exp3times0_25vintchl,
':',color='violet',lw=3,label='* 1/4')
h0_50=ax.plot(first_year,array1d_iT_exp3times0_50vintchl,
'-',color='blue',lw=3,label='* 1/2')
h1_00=ax.plot(first_year,array1d_iT_modvintchl,
'-',color='black',lw=3,label='default')
h2_00=ax.plot(first_year,array1d_iT_exp3times2_00vintchl,
'-',color='orange',lw=3,label='* 2')
h4_00=ax.plot(first_year,array1d_iT_exp3times4_00vintchl,
':',color='red',lw=3,label='${[NO_3]}^{winter}$ * 4')
ax.set_yscale('log')
ymin=1E0
ymax=1E2
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('Vertically integrated Chl $a$\n($\mathrm{ mg\ Chl\ m^{-2} }$)')
ax.yaxis.set_ticks_position('both')
ax.set_ylim(ymin,ymax)
ax.grid()
#get handles and labels
handles, legend_labels = ax.get_legend_handles_labels()
#specify order of items in legend
order = [2,4,3,1,0]
#add legend to plot
ax.legend([handles[idx] for idx in order],
[legend_labels [idx] for idx in order],
loc='upper left',
bbox_to_anchor=(0.08,1.00))
ax.axvline(obs_breakup,color='k',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
plt.text(0.03,0.90,'(b)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- EXP-0: ADDITIONAL X-AXIS
ax=axs[2]
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels(labels)
ax.axes.get_yaxis().set_visible(False)
ax.spines['top'].set_visible(False)
# --- EXP-0: SCHEMATIC
ax=axs[3]
h0=ax.plot(np.array(df.exp0x),np.array(df.exp0y),
'-',color='black',lw=3,label='EXP-0')
ax.set_xlabel('Time')
ax.xaxis.set_label_coords(0.05,-0.05)
ax.set_xlim(xlims)
xticks=([df.exp0x.iloc[2]])
xticklabels=[r'$\leftarrow t_P \rightarrow$']
ax.set_xticks(xticks)
ax.set_xticklabels(xticklabels)
ax.set_yscale('log')
ax.set_ylabel('Log10(Chlorophyll $a$)')
ax.yaxis.set_label_coords(-0.075,0.50)
ymin=1E-4
ymax=300
ax.set_ylim(ymin,ymax)
ax.tick_params(axis='y',which='minor',left=False,labelleft=False)
yticks=([df.exp0y.iloc[2]])
yticklabels=[r'$\uparrow$'
'\n'
r'$y_P$'
'\n'
r'$\downarrow$']
ax.set_yticks(yticks)
ax.set_yticklabels(yticklabels)
ax.plot(df.exp0x.iloc[1],df.exp0y.iloc[1],
'o',color='black')
ax.text(df.exp0x.iloc[1]-6,df.exp0y.iloc[1]+1,'I')
ax.plot(df.exp0x.iloc[2],df.exp0y.iloc[2],
'o',color='black')
ax.text(df.exp0x.iloc[2]-6,df.exp0y.iloc[2]+20,'P')
ax.axvline(df.exp0x.iloc[2],
ymax=( math.log10(df.exp0y.iloc[2])-math.log10(ymin) ) \
/ ( math.log10(ymax) -math.log10(ymin) ),
color='black',linestyle='--')
ax.axhline(df.exp0y.iloc[2],
xmax=df.exp0x.iloc[2]/xmax,
color='black',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
ax.legend(loc='upper right')
plt.text(0.03,0.85,'(c)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
plt.text(0.5,1.05,'Standard magnitude of the bloom peak',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- EXP-3: SCHEMATIC
ax=axs[4]
h0=ax.plot(np.array(df.exp3x),np.array(df.exp3y),
':',color='violet',lw=3,label='EXP-3')
ax.set_xlabel('Time')
ax.xaxis.set_label_coords(0.05,-0.05)
ax.set_xlim(xlims)
xticks=([df.exp3x.iloc[2]])
xticklabels=[r'$\leftarrow t_P \rightarrow$']
ax.set_xticks(xticks)
ax.set_xticklabels(xticklabels)
ax.set_yscale('log')
ax.set_ylabel('Log10(Chlorophyll $a$)')
ax.yaxis.set_label_coords(-0.075,0.50)
ymin=1E-4
ymax=300
ax.set_ylim(ymin,ymax)
ax.tick_params(axis='y',which='minor',left=False,labelleft=False)
yticks=([df.exp3y.iloc[2]])
yticklabels=[r'$\uparrow$'
'\n'
r'$y_P$'
'\n'
r'$\downarrow$']
ax.set_yticks(yticks)
ax.set_yticklabels(yticklabels)
ax.plot(df.exp3x.iloc[1],df.exp3y.iloc[1],
'o',color='black')
ax.text(df.exp3x.iloc[1]-6,df.exp3y.iloc[1]+1,
'I\'\'\'',color='violet')
ax.plot(df.exp3x.iloc[2],df.exp3y.iloc[2],
'o',color='black')
ax.text(df.exp3x.iloc[2]-6,df.exp3y.iloc[2]+20,
'P\'\'\'',color='violet')
ax.axvline(df.exp3x.iloc[2],
ymax=( math.log10(df.exp3y.iloc[2])-math.log10(ymin) ) \
/ ( math.log10(ymax) -math.log10(ymin) ),
color='black',linestyle='--')
ax.axhline(df.exp3y.iloc[2],
xmax=df.exp3x.iloc[2]/xmax,
color='black',linestyle='--')
rect=mpl.patches.Rectangle(xy=(xmin,ymin),
width=obs_breakup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
rect=mpl.patches.Rectangle(xy=(sim_freezeup,ymin),
width=364-sim_freezeup,height=ymax,
linewidth=0,facecolor='grey',alpha=0.25)
ax.add_patch(rect)
ax.legend(loc='upper right')
plt.text(0.03,0.85,'(d)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
plt.text(0.5,1.05,'Low magnitude of the bloom peak',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- POSITION
axs[0].set_position([0.13,0.92, 0.85,0.03])
axs[1].set_position([0.13,0.66, 0.85,0.24])
axs[2].set_position([0.13,0.64, 0.85,0.01])
axs[3].set_position([0.13,0.34, 0.85,0.24])
axs[4].set_position([0.13,0.04, 0.85,0.24])
with plt.style.context('mplstyles/exp.mplstyle'):
# Plot
fig,axs=plt.subplots(5,1,figsize=(16,19.5))
make_plots(axs)
# --- SAVE
plt.savefig('figures_progress/exp3.png',dpi=1000)
The conceptual model and the schematic.
def make_plots(axs):
xlims=(first_year[0],first_year[-1])
xmin=xlims[0]
xmax=xlims[1]
# --- CONCEPTUAL MODEL
ax=axs[0]
im=PIL.Image.open('figures_progress/conceptual_model.png')
ax.imshow(im)
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)
plt.text(0.03,0.86,'(a)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- SCHEMATIC
ax=axs[1]
h0=ax.plot(np.array(df.exp0x),np.array(df.exp0y),
'-',color='black',lw=3,label='EXP-0')
ax.set_xlabel('Time')
ax.xaxis.set_label_coords(0.05,-0.03)
ax.set_xlim(xlims)
xticks=([df.exp0x.iloc[1],df.exp0x.iloc[2]])
xticklabels=[r'$\leftarrow t_I \rightarrow$',
r'$\leftarrow t_P \rightarrow$']
ax.set_xticks(xticks)
ax.set_xticklabels(xticklabels)
ax.set_yscale('log')
ax.set_ylabel('Log10(Chlorophyll $a$)')
ax.yaxis.set_label_coords(-0.075,0.50)
ymin=1E-4
ymax=300
ax.set_ylim(ymin,ymax)
ax.tick_params(axis='y',which='minor',left=False,labelleft=False)
yticks=([df.exp0y.iloc[0],df.exp0y.iloc[2]])
yticklabels=[r'$\uparrow$'
'\n'
r'$y_I$'
'\n'
r'$\downarrow$',
r'$\uparrow$'
'\n'
r'$y_P$'
'\n'
r'$\downarrow$']
ax.set_yticks(yticks)
ax.set_yticklabels(yticklabels)
ax.plot(df.exp0x.iloc[1],df.exp0y.iloc[1],
'o',color='black')
ax.text(df.exp0x.iloc[1]-6,df.exp0y.iloc[1]+1,'I')
ax.plot(df.exp0x.iloc[2],df.exp0y.iloc[2],
'o',color='black')
ax.text(df.exp0x.iloc[2]-6,df.exp0y.iloc[2]+20,'P')
ax.axvline(df.exp0x.iloc[1],
ymax=( math.log10(df.exp0y.iloc[1])-math.log10(ymin) ) \
/ ( math.log10(ymax) -math.log10(ymin) ),
color='black',linestyle='--')
ax.axvline(df.exp0x.iloc[2],
ymax=( math.log10(df.exp0y.iloc[2])-math.log10(ymin) ) \
/ ( math.log10(ymax) -math.log10(ymin) ),
color='black',linestyle='--')
ax.axvline(df.exp0x.iloc[2],
ymax=( math.log10(df.exp0y.iloc[2])-math.log10(ymin) ) \
/ ( math.log10(ymax) -math.log10(ymin) ),
color='black',linestyle='--')
ax.axhline(df.exp0y.iloc[2],
xmax=df.exp0x.iloc[2]/xmax,
color='black',linestyle='--')
ax.axhline(df.exp0y.iloc[0],
xmin=df.exp0x.iloc[1]/xmax,
xmax=df.exp0x.iloc[2]/xmax,
color='black',linestyle='--')
# ref.:
# https://stackoverflow.com/questions/68111673/matplotlib-draw-ellipse-to-annotate-on-top-of-logarithmic-scale
# answer by tmdavison
# Arc centre coordinates
x,y=df.exp0x.iloc[1],df.exp0y.iloc[1]
# use the axis scale tform to figure out how far to translate
arc_offset=mpl.transforms.ScaledTranslation(x,y,ax.transScale)
# construct the composite tform
arc_tform = arc_offset + ax.transLimits + ax.transAxes
# Create the arc centred on the origin, apply the composite tform
width,height=20,1 # not exact, by trial and error
theta1,theta2=0,3 # not exact, by trial and error
arc=mpl.patches.Arc(xy=(0,0),width=width,height=height,
color="grey",fill=False,lw=1,
theta1=theta1,theta2=theta2,
transform=arc_tform)
ax.add_patch(arc)
ax.text(174,1.2,'$m$')
plt.text(0.03,0.86,'(b)',
horizontalalignment = 'left',verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- POSITION
axs[0].set_position([0.13,0.50, 0.85,0.50])
axs[1].set_position([0.13,0.04, 0.85,0.46])
with plt.style.context('mplstyles/exp.mplstyle'):
# Plot
fig,axs=plt.subplots(2,1,figsize=(16,10.2))
make_plots(axs)
# --- SAVE
plt.savefig('figures_progress/discussion.png',dpi=1000)